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Abstract 


In order to calculate the velocity and passive scalar fields in 
homogeneous turbulent flows , the three-dimensional , time-dependent 
equations of motion and the diffusion equation have been solved numer- 
ically. The following cases have been treated: 

1. isotropic, homogeneous turbulence with decay of a passive scalar, 

2. homogeneous turbulent shear flow with a passive scalar whose mean 

varies linearly in the spanwise direction. 

The solutions have been obtained at relatively low Reynolds numbers so 
that all of the turbulent scales could be resolved without modeling. 

The numerical simulations were carried out on the ILLIAC-IV. The 
computeir program used 16 x 16 x 16, 32 x 32 x 32, or 64 x 64 x 64 
uniformly spaced mesh points. The pseudo-spectral method was used for 
space differencing and the second-order Adams-Bashforth method for time 
differencing. 

Turbulent statistics such as integral length scales , Taylor micro- 
scales, Kolmogorov length scale, one- and two-point correlations of 
velocity-velocity and velocity-scalar, turbulent Prandtl/Schmidt number, 
r.m.s. values of velocities, the scalar quantity and pressure, skewness, 
decay rates, and decay exponents have been calculated. Also, one- and 
three-dimensional spectra of velocity components and scalar quantity, 
three-dimensional spectra of production, dissipation and transfer terms, 
pressure, Reynolds stress, and heat flux terms have been calculated. 
Some of the results are compared with the available experimental re- 
sults, and good agreement is obtained. 

The behavior of the Reynolds stress And the turbulent heat flux 
have been carefully studied. The conventional models for the turbulent 
Reynolds stress and heat flux, turbulent Prandtl number, and pressure- 
strain and pressure-scalar covariances have been tested, and constants 
for the models have been obtained. 
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Chapter I 


INTRODUCTION 


1 . 1 General Background 

For a century, turbulence has been the subject of intensive re- 
search, both theoretically and experimentally. Due to its complexity 
and despite strong efforts , turbulence remains one of the most chal- 
lenging unsolved problems of the physical sciences. History suggests 
that the road to the general solution of the problem is very long , if 
not endless. 

Turbulence originates in instabilities of laminar flow. Generally, 
the instabilities produce wave-like structures which can absorb energy 
from the mean flow. As they grow, nonlinear effects cause energy trans- 
fer to other modes , and eventually the noisy patterns that are generally 
regarded as "turbulence" results. Energy obtained from the mean flow is 
eventually transferred to smaller-scale structures , where it will be 
dissipated. 

Turbulence contains structures with a wide range of scales. The 
large-scale turbulent structures differ markedly from one flow type to 
another, while the small scales tend to be more universal. Since nei- 
ther experimentalists nor computors have been able to resolve all of the 
turbulent structures for high Reynolds number flows, approximations have 
had to be made. 

Almost all approaches to the prediction of turbulent flows are 
based on Osborne Reynolds' idea of averaging the Navier-Stokes equations 
over an ensemble of identical flows, or some equivalent (time or span 
averaging) to obtain equations for the mean velocities. Due to the 
nonlinearity of the Navier-Stokes equations, the averaged equations are 
not closed. Therefore, turbulence models are needed to bridge the gap 
between the difficulty of the equations and our limited computational 
capability to handle the equations. 

Turbulence modeling has made considerable strides in the past two 
decades. The testing of the models has, until recently, been a highly 
empirical affair. It is currently impossible to measure all of the 
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factors needed to test the models experimentally. However, it is pos- 
sible to use exact numerical solutions of the equations of motion for 
turbulent flows as the basis for testing models. The difficulty here is 
that , with the available computers , one can compute only relatively low 
Re3molds number flows. In addition, only a few simple turbulent flows 
can be computed with sufficient accuracy. Due to the large range of 
scales of turbulent structures, we cannot expect to compute high Rey- 
nolds number turbulent flows directly with the next few generations of 
computers. Despite these limitations, one might compute relatively 
simple flows at moderate Reynolds numbers and extrapolate the results to 
higher Reynolds numbers. 

The simplest type of turbulence is homogeneous isotropic turbu- 
lence. This is defined as a flow whose statistical features have no 
preference for any specific direction or location. Thus a minimum 
number of quantities and relations are required to describe its struc- 
ture and behavior. Although it is an Idealized flow, its characteris- 
tics form a basis for the study of more complex non-isotroplc turbulent 
flows. Moreover, in actual turbulent flows at high Reynolds numbers, 
the fine structures seem to be nearly isotropic, so the results of 
isotropic turbulence may be useful for characterizing the small scales 
of turbulence. Because of its simplicity and importance, extensive 
theoretical work has been done on this subject and a number of exper- 
iments have been carried out. 

A somewhat more complex type of turbulent flow is homogeneous shear 
flow. The role of shear in most turbulent flows is to couple the mean 
velocity field with the turbulent velocity field. Homogeneous shear 
flow is a problem of complexity intermediate between unsheared homogene- 
ous and Inhomogeneous turbulent shear flows. It has been generated both 
experimentally and computationally. Its study can provide better under- 
standing of turbulence and can be used in testing the validity of turbu- 
lence theories and statistical turbulence models. 

It is useful to consider the case in which a scalar quantity is 
present in a turbulent flow. The scalar may be temperature or en- 
thalpy, or, in a mixture of fluids, species concentration. The study of 
the scalar field has attracted less attention than the study of the 
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hydrodynamics itself. However, the investigation of turbulent mixing is 
£in area of great importance. It is required in any situation that in- 
cludes mass or heat transfer, for example, in the atmosphere. Turbulent 
mixing is also important in reacting flows and combustion. Prediction 
of most turbulent reacting and combusting flows is essentially impos- 
sible without accurate treatment of mixing. Study of a passive scalar 
in turbulent flows will help us understand the turbulent mixing proc- 
esses that occur in these flows. 

One of the objectives of the present report is to compute three- 
dimensional, time-dependent solutions for isotropic turbulence and 
homogeneous shear turbulent flows by numerical methods, and to study the 
scalar field as well as the velocity field. The results will be used in 
a number of ways . 

1.2 Survey of Literature 

The mathematical development of hydrodynamics was started in the 
eighteenth century. The dynamic equations for inviscid flows were 
formulated mainly by Euler (1752). The equations of motion for a vis- 
cous flow were derived by Navier (1827), Poisson (1831), Saint-Venant 
(1843), and Stokes (1845). 0. Reynolds (1883) introduced the non- 
dimensional Reynolds number and used it as a criterion to distinguish 
laminar flows from turbulent flows. Boussinesq (1877) introduced the 
concept of eddy viscosity. Mixing-length theory was introduced by 
Prandtl (1925) and von Karman (1930). 

A. Homogeneous Isotropic Turbulence 

The concept of isotropic turbulence was Introduced by Taylor 
(1935). Von Karman (1937) introduced the use of tensors in isotropic 
turbulence, while von Karman and ibwarth (1938) studied the statistical 
theory of turbulence. In 1941, Kolmogoroff suggested that the small- 
scale components of turbulence are approximately isotropic, while Batch- 
elor and Stewart (1949) showed that large scales of the turbulence gen- 
erated behind a grid are anisotropic. 

Almost all experiments on isotropic turbulence use rigid, uniform 
grids to generate the turbulence. The first successful attempt to 
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generate nearly isotropic turbulence was that of Simmon and Salter 
(1934). They found that at high Reynolds numbers the turbulence far 
behind a grid is a good laboratory realization of isotropic turbu- 
lence. Several other experiments have been carried out for this case, 
among which that of Comte-Bellot and Corrsin (1971) provides rather 
complete statistical information. Frenkiel, Klebanoff, and Huang (1979) 
carried out the experiment in both water and air. 

The final stage of isotropic turbulence was studied experimentally 
by Batchelor and Townsend (1950), Tan and Ling (1963), Lee (1965), and 
Bennett and Corrsin (1978). Also Tavoularis, Bennett, and Corrsin 
(1979) studied the velocity skewness of isotropic turbulence at small 
Reynolds numbers . 

Isotropic homogeneous turbulence has been studied by computers. 
Numerical simulation of the decay of two-dimensional, isotropic, 
homogeneous turbulent flows was carried out by Herring, Orszag, Kraich- 
nan, and Fox (1974), and turbulence theories were tested. The first 
simulation of three-dimensional, homogeneous, isotropic turbulence was 
made by Orszag and Patterson (1971). Schumann (1975) also simulated 
isotropic turbulence numerically and studied the velocity and pressure 
fluctuations. Three-dimensional, homogeneous, isotropic. Incompressible 
turbulence has been studied at Stanford University using the large eddy 
simulation technique. Kwak, Reynolds, and Ferziger (1975), and Shaanan, 
Ferziger, and Reynolds (1975) simulated isotropic turbulence and ob- 
tained very good agreement with the experimental results of Comte-Bellot 
and Corrsin. Clark, Ferziger, and Reynolds (1977) carried out a direct 
calculation for the isotropic, homogeneous turbulent flows at low Rey- 
nolds number and used it to study subgrid scale modeling. 

There is also a series of experiments in heated isotropic turbu- 
lence. In all cases the heating was small enough that the temperature 
could be assumed a passive scalar and density variations neglected. The 
diffusion of heat from a fixed line source in grid-generated , nearly 
isotropic turbulence was studied experimentally by Schubauer (1935) , 
Collis (1948), Frenkiel (1950), Townsend (1951), Uberoi and Corrsin 
(1952), Shlien and Corrsin (1974), and Libby (1975). These experiments 
studied the downstream development of temperature fluctuations and 
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measured statistical properties such as length scales, decay rates, 
velocity-velocity and velocity-temperature correlations, and spectra. 

The fluctuating temperature field in isotropic turbulence generated 
by uniformly heated grids was studied by Kistler, O'Brien, and Corrsin 
(1956), Mills and Corrsin (1959), Yeh and van Atta (1973), Sepri (1976), 
and Warhaft and Lumley (1979). They have presented various correla- 
tions, length scales, and spectra of the scalar and velocity fields. 
Table 1-1 summarizes the basic parameters of the above experiments. Yeh 
et al. (1974) carried out their experiment at relatively low Reynolds 
number and studied the spectral transfer of the scalar and velocity 
fields. Table 1-1 shows that, unlike the decay exponent of the turbu- 
lent kinetic energy, the decay exponent of the fluctuating scalar inten- 
sity varies considerably from one experiment to another. Domls (1979) 
simulated the flow by large eddy simulation and studied the effects of 
the initial length scales on the decay exponent of the scalar. Ife con- 
cluded that the decay exponent of the scalar is a linear function of the 
Initial length scale ratio. We shall study the effects of Reynolds 
number and Prandtl number on the decay exponents. 

v/e« 

Non-unlformly heated grids were used by Wiskind (1962) and Venka - l ' t if- 
and Chevray (1978) to generate a uniform temperature gradient in 
isotropic turbulence. Measurements of probability and joint probability 
densities of velocity and temperature as well as spectra and the other 
turbulence statistics were obtained in the latter experiment. Cornelius 
and Foss (1978) measured the diffusion of particles in isotropic turbu- 
lence by a unique method. 

B. Homogeneous Shear Flows 

Reis (1952) studied homogeneous shear flows mathematically. Rose 
(1966) experimentally generated this flow and studied the behavior of 
the turbulence Intensities. Champagne, Harris, and Corrsin (CHG) 
(1970), Harris, Graham, and Corrsin (HGC) (1976), Rose (1970), Ife^ang 
(1971), and Mulhearn and Luxton (1975) studied this flow at higher 
Reynolds numbers. These experiments measured the shear stress correla- 
tions , various length scales , one- and two-point velocity correlations , 
spectra, and the kinetic energy of turbulence. Table 1-2 summarizes 
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some of the basic parameters of these flows. Rose's experiment was 
carried out at lower Reynolds number and the Harris et al. experiment at 
higher shear rate. 

Numerical simulations of homogeneous shear flows were made by 
Shannan, Ferzlger, and Reynolds (1975) using large eddy simulation. 
Later, R. S. Rogallo (1977) made direct simulations of Reynolds number 
homogeneous turbulent shear flows. Similar calculations were done for 
the compressible case by Felerelsen, Reynolds, and Ferzlger (1981). 

Temperature as a passive scalar in homogeneous shear flows was 
studied by Tavoularis (1978). Hte generated a turbulent homogeneous 
shear flow with a uniform mean temperature gradient and studied the 
statistics of the scalar field, as well as the velocity field. The 

basic parameters of this experiment are similar to those of the Harris 
et al. experiment and are shown in Table 1-2. He used air (Prandtl 
number 0.7) and the mean- temperature gradient was 9.5 C/m. We shall 
compare our computed scalar field results with this experiment. Fox 
(1964) simulated the heat transfer in a two-dimensional, homogeneous 
shear flow numerically, and studied the effect of molecular Prandtl 
number on the turbulent heat flux and the turbulent Prandtl number. 

1.3 ^fc>tlvations and Objectives 

Homogeneous shear flow is a flow of complexity between that of 
unsheared isotropic homogeneous turbulence and inhomogeneous shear 
flows. Although simple (it does not contain rigid boundaries, laminar- 
turbulent interfaces, or other inhomogeneities), it retains many basic 
features of sheared turbulence, such as production of turbulent kinetic 
energy from mean shear and non-zero Reynolds shear stress. This flow 
has been generated experimentally by several investigators mentioned in 
the previous section. From this flow, many turbulent statistics can be 
calculated and turbulence theories and models tested. 

The investigation of mixing of a scalar in turbulent flows is of 
great Importance. Turbulent mixing occurs in any flow which includes 
mass or heat transfer. There are many applications in meteorology and 
oceanography, as well as engineering. The last includes reacting flows 
as well as flows in which heat transfer is important. 


6 




We shall study mixing of a passive scalar in homogeneous turbulent 
shear flows and compare the scalar field with the velocity field. This 
helps us to understand some basic phenomena of turbulence. By means of 
analogy between the velocity field and the scalar field, we shall study 
various aspects of turbulent mixing. We shall also test some turbulence 
models . 

Our specific objectives are these: 

1. To compute the mixing of a passive scalar in isotropic turbulent 
flow at relatively low Reynolds number by solving the Navler-Stokes 
and concentration equations numerically. We shall calculate r.m.s. 
values of velocity components and the scalar quantity, various 
length scales, spectra, and the decay exponent ratio between the 
scalar field and the velocity field. We shall study the effects of 
the Reynolds and Prandtl numbers on the decay rate of turbulent 
kinetic energy and the r.m.s. value of the scalar quantity. All of 
these results will be compared with experiment. 

2. To compute the behavior of a passive scalar in a turbulent homo- 
geneous shear flow, with uniform scalar gradient. In this flow we 
shall calculate one- and two-point correlations, turbulent Prandtl/ 
Schmidt number, r.m.s. values of velocity, scalar, and pressure, 
skewness, various length scales, and spectra. We shall study the 
effects of Reynolds number, Prandtl number, and shear number (a 
non-dimensional shear rate) on the turbulence statistics. We shall 
present correlations for many of the turbulence statistics. Spe- 
cial attention will be given to the turbulent Reynolds stresses and 
the turbulent heat flux. Finally, energy and heat transfer between 
the various sizes of eddies will be discussed by calculating the 
spectra of the transfer terms. 

3. To test various models for heat/ mass flux, Reynolds stress, 
pressure-strain and pressure-scalar covariences. The validity of 
the models as well as the constants for each models will be 
studied. 



Table 1-1 


Descriptions and Results of Various Experiments and Numerical 
Simulations on Decay of Temperature Fluctuations 
in Grid-Generated Homogeneous Isotropic Flows 



Rex 

„ ** 

Pr 

n 

m 

^u^^e 

Yeh & Van 
Atta (1973) 

60.97 

10550 

0.725 

1.37 

1.33 

0.718 

Sepri (1976) 

60.62 

10484 

0.72 

1.39 

1.4 

? 

War haft and 
Lumley (1978) 

130 

10000 

0.72 

1.34 

.858-1.407 

.216-0.316 

Domis (1979) 

80.97 

10550 

0.72 

0.885 

0.91-2.03 

0.71-0.93 


Rei^ is the Reynolds number based on the Taylor microscale and turbu- 
lent velocity. 


A 

Re^^ is the Reynolds number based on the mesh size and streamwise 
velocity. 
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Table 1-2 


Description of Various Experiments 
on Homogeneous Shear Flows 



Shear Rate 

Reynolds No. 

Shear No. 

Total Shear 

Experiment 

S (sec”^) 

r, * 

Re^ 


St 

Rose (1966) 

13.68 

133.0 

9.53 

3.2 

Champagne et al. 

12.90 

204.8 

7.79 

3.2 

(1970) 





Harris et al. 

44.0 

396.4 

8.28 

10.81 

(1977) 





Tavoularis (1979) 

46.8 

397.9 

14.76 

12.65 


The values of Re^^ and ^ are given at Xj^/h = 7.5 for Rose and 
Harris et al., at Xj^/h = 8.5 for Champagne et al., and at Xj^/h = 1.1 
for Tavoularis. 

The shear number is defined , as = SL/q , where L is integral 

length scale and q = < 
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Chapter II 


COMPUTATIONAL METHODS 


One of the major difficulties in numerical calculations of turbu- 
lent flows is the wide range of length scales in the flow. To capture 
the large scales of turbulence , the computational box has to be bigger 
than the largest scale. On the other hand, to resolve the smallest 
scale of turbulence in a grid-based calculation, the mesh size has to be 
smaller than the smallest length scale, which is the Kolmogoroff micro- 
scale defined as: 


n 



( 2 - 1 ) 


where v is the kinematic viscosity and e is the dissipation rate, 
i.e., the rate of destruction of turbulent kinetic energy by viscosity. 
The dissipation is also related to the large scale by the relation: 

e = C^q^/Ld (2-2) 

where q is the r.m.s. velocity (q^ = < u^ > + < u| > + < u| > ) , 
is a large eddy length scale, and is a proportionality constant. 

One can find the number of mesh points required to resolve both the 
large and small scales in a three-dimensional computation by combining 
the above relations. The number of mesh points required for accurate 
resolution is approximately 



and depends strongly on Reynolds number. Since, from a practical point 
of view, the number of mesh points is limited by the computational speed 
and size of the computer memory, computations are currently limited to 

O 

low turbulent Reynolds numbers. In our case, with 64"^ cubed mesh 
points, the Reynolds number has to be less than about 100. 
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2.1 Governing Equations 


The equations of motion for an Incompressible fluid are: 


3U. 


X 


3T 


3C 

3T 


3U.U . 
i J 

3P 

a- 

(2-4) 

3X, 

3X^ 

3X.3X, 

j 

3U. 

1 

J 3 


^ ■ 

3CU. 

1 

3Xi 

0 


(2-5) 

+ D 

a C 

3X.X . 
J 3 

(2-6) 


where U, P, v, C, and D are velocity, ratio of pressure to density, 
kinematic viscosity, scalar concentration, and diffusion coefficient, 
respectively. The summation convention Is Implied. Note that the sca- 
lar quantity can be temperature or concentration of a species. 


Instead of satisfying the Incompressible continuity equation, we 
use the Poisson equation. This equation Is derived by taking the diver- 
gence of Eq. (2-4) and using Eq . (2-5) to eliminate some of the terms. 
The result Is 

2 

'''■ - --SX-M (2-7) 

i J 


2 .2 Numerical Approximations 

The governing equations given above do not have analytical solu- 
tions for turbulent flows, and numerical approximations must be made. 
Since the data will be discretized onto mesh points , we must approximate 
partial derivatives In terms of the data at these points. We use the 
pseudo-spectral method for this purpose. This method, which Is good for 
problems with periodic boundary conditions, has been analyzed by Lanczos 
(1956) and has been proposed by Orszag and Patterson (1971) and Krelss 
and Ollger (1973) as an approximation method and used by many computers 
since then. Mansour et al. (1977) used the pseudo-spectral method and 
compared It with some second-order and fourth-order finite-difference 
schemes. They showed that use of the pseudo-spectral method Improved 
accuracy considerably. The method is briefly explained below. 


11 



Suppose a function of one variable F(X) is periodic with period 
L. For equally spaced mesh points such a function can be approximated 
at the mesh points as follows: 


. ik X. 

F(X,) = 2^ F(k^) e J 

n=-N/2 




where 


X 


N 

3 

h 


k = 


number of mesh points , 

h;j , j = 0, 1,2,3, ... ,N-1, 

mesh size, 

2 IT 


n Nh 


- N N , , 

n, n u,...,“ i, and 




-ik X . 

a *«• 


(2-8) 


(2-9) 


To compute F(k) , we use the Fast Fourier Transform (FFT) , origin- 
ally developed by Cooley and Tukey (1963). This algorithm is good for N 
= 2^^, where n is an integer (see, for example, Bracewell (1965). 
Newer methods are efficient for other numbers of points but are not 
needed for our purposes. 


If we regard the expansion (2-8) as an interpolation formula and 
replace Xj by a continuous variable X, the derivative of F(X) with 
respect to X can be approximated by 

n 

Therefore, to compute dF(X)/dX, we calculate the Fourier transform 
of F(X) , multiply it by ik, and take its inverse Fourier transform. 

To advance in time, a second-order Adams-Bashforth (AB) method was 
used. It has been shown by Billy (1965) that this method is weakly 
unstable, but its spurious computational production of kinetic energy is 
small. A fourth-order Runge-Kutta method would be more stable; however, 
it is more expensive and requires larger computer memory. 

The AB method of advancing the solution of 9U/3t = H from the 
time step n to time step nH-1 is 
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( 2 - 11 ) 


yii+l ^ i At(3H° - 

where At Is the time step. 

As we see, it is a two-step explicit method and needs a starting 
method; we use the Euler method to start the calculation. The Euler 
formula is 


= U’^ + AtH’^ (2-12) 

2.3 Initial Conditions 

We start with an initial velocity field that is randomly chosen, 
but which satisfies continuity and has a given three-dimensional spec- 
trum. To obtain such an Initial field, we first use computer-generated 
random numbers for each component of velocity at each mesh point, . 

To this we add the gradient of a potential which will be chosen to make 
the field divergence-free; 





(2-13) 


To make 

equation and set 


divergence-free, we take the divergence of the above 
to find 


= 0 


a ^ 
Sx^ax^ 


au/^> 

W. 

1 


(2-14) 


This is a Poisson equation for the potential and is solved by a Fourier 
transform method. The Fourier transform of Eq. (2-14) is: 


1 '' 

k <j) = 


- 


(2-15) 


From this equation, we calculate <(), Invert the Fourier transform to 
get (|), and substitute into Eq. (2-13) to get the desired divergence- 
free field. 

Now we take the Fourier transform of the resulting field and mul- 
tiply each mode by a factor such that the desired three-dimensional 
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spectrum is obtained. The resulting field is nearly Isotropic, 
divergence-free, and has the desired spectrum. 

We use a similar procedure to calculate the initial scalar field. 
We start with a field of random numbers and take its Fourier transform. 
Then we multiply each mode by a factor which yields the desired spectrum 
for the passive scalar. 

The Initial three-dimensional spectra of the velocity and scalar 
fields are step functions. The width, location, and height of the spec- 
tra determine the initial Taylor microscales, integral length scale, and 
the r.m.s. values of turbulent velocity components and scalar concentra- 
tion. All the initial values are non-dimensionalized (the computational 
box length is 2u) and must be scaled to obtain the dimensional values. 
The problem contains three dimensional units: length, time, and tempera- 
ture (or mass in the case where the scalar is concentration). To dlmen- 
sionalize the results, we need a length, a time, and a temperature or 
mass, or three independent combinations of these parameters. 

2 . 4 C omputation of Ebmogeneous Isotropic Turbulent Flows with Decay of 

a Passive Scalar 

In this section we consider a method of computation for simulating 
a grid-generated isotropic homogeneous flow. The computational box is 
shown in Fig. 2-1. 

The governing equations for this case are Eqs. (2-4), (2-6), and 
(2-7). The mean values of velocity, pressure, and scalar quantity are 
set to zero. The Initial spectra of velocity and the scalar are given 
in Fig. 2-2. For all cases we use 32 x 32 x 32 equally spaced mesh 
points. The accuracy and stability of the numerical method are dis- 
cussed in Appendix A and will be explained later in Section 3D of this 
chapter. 

The computing time per time step is approximately 3.2 seconds of 
CPU time on the ILLIAC IV. 
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2.5 Computation of Homogeneous I\irbulent Shear Flow with Decay of a 
Passive Scalar 


In this section we derive governing equations numerically to simu- 
late turbulent homogeneous shear flow with a passive scalar. We also 
describe the computational details. 


A. Mean Values of Velocities and Scalar Quantity 

In homogeneous flows , the mean values of the velocity components 
and scalar concentration may vary with time and/or space. We are inter- 
ested In mixing a passive scalar in a homogeneous shear flow with the 
following mean values of velocity, pressure, and concentration: 

Ui = SX2 

?2 = 0 

?3 = 0 (2-16) 

P = 0 

^ “ ^ 2^2 

where an overbar denotes an ensemble average and S and S 2 are shear 
rate and mean concentration gradient, respectively. The mean values, as 
shown above, remain unchanged throughout the simulation. 

The fluctuating quantities are defined as follows: 


or 


and 



(2-17) 


(2-18) 
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B. Coordinate Transformations 


If we substitute Eqs. (2-15) and (2-16) into Eqs. (2-4) 
(2-6), we get the set of equations 

9u, 


(2-5) , and 


9X 


= 0 


9u^ 

9t 


9u .u . 
1 J 


i 

9u, 


9X, 

J 




3 u. 


+ V 


96 

9t 


96u, 


9X.9X. 
J J 


9X 


- SX, 


96 ^ ^ „ 3^6 

2 3X, ^2^^2 ° 9X.9X . 

1 J J 


(2-19) 


Periodic boundary conditions cannot be applied to these equations. 
They are prohibited by the presence of SX 2 as a coefficient in one 
term of each of the last two equations. These terms are different at 
the top and bottom of the computational box, and application of periodic 
boundary conditions to them would result in numerical errors similar to 
the Gibbs phenomenon at these two boundaries. We therefore introduce a 
coordinate transformation which eliminates the non-periodic terms. The 
new coordinates are 


Xj^ = - SX2t 

X 2 = X 2 (2-20) 

The new coordinate system moves with the mean flow, as shown in Fig. 
2-3. 

Transforming the equations to the new coordinate is a little com- 
plicated and is done by using the chain rule. Details are given in 
Appendix B. The resulting equations will be given below. 

C . The Governing Equations 

The governing equations are obtained by applying the transformation 
(2-20) to equations (2-19). The transformed continuity equation is 

Ujj^ j - St U 2 = 0 (2-21) 
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The momentum equations become 


9u. 


- = - (u.u.) . + St u.u„ - Su„6 - p + v(u. ..-2Stu. , ,.+S^t^u . 


JT 


i j ,j 


1 2 


2 il 


1 ,JJ 


'i,12' 


1,11^ 


(2-22) 


The Poisson equation for the pressure is 


P .. - 2St p St'^p = H, . - St a, , 
,ix *^,12 *^,11 i ,1 z,l 


(2-23) 


Finally, the scalar equation becomes 


If ■ - - ‘*2"2 + " I 


'.u_ 


In Eq. (2-22) 


\ = ~ ■ ^“2^i2 ^i.l2+^^*^S,ll^ 

(2-25) 


The governing equations for turbulent homogeneous shear flow with 
mixing of a passive scalar are (2-22) through (2-24) . The transformed 
equations admit periodic boundary conditions, which in turn permit the 
use of spectral methods of computing derivatives . 


The Poisson equation for the pressure must be derived such that the 
continuity equation is satisfied at the following time. The discretized 
equation for pressure which has this property when we use the Adams- 
Bashforth method for the time differencing has the following form.: 


3xj 3xj " S(2n+1) At 


2 2 
3 ^ , ^,,.2 3 

+ S n( n+1 ) A — ; 


3x! 


1 -1 


P - “i,i 


- S(n+l)AtH2 1 - J 


-.n _ 


- p^u p^u 


- S(n+1) 


P n-1 n-1 

[_^,1 " P,12 


+ S(n-l) Atp At 


n-1-] 


(2-26) 


(2-27) 
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D . Remeshlng 


There are terms in Eqs. (2-21) through (2-24) which contain St as 
a coefficient. As time increases, these terms become large and cause 
numerical problems. In other words, the and x -2 axes become more 

nearly parallel as time increases and the equations become stiff. To 
avoid this problem we use a procedure which was first used by Rogallo 
(1978). We start the calculation with coordinates orthogonal (Fig. 2- 
3a). Then, after a number of time steps (see Fig. 2-3b) we stop the 
calculation and rotate the Xj^ = 0 line to the other side of the 
normal, as shown in Fig. 2-3c, and continue the calculations. This 
procedure is called remeshing and is exact only if it is done at St = 
0.5, 1.0, or 1.5, etc. for these times, there is no need to interpolate 
the data to get the information at the new mesh points after 
remeshlng . However , there is an aliasing problem associated with 
this. To show this, let us assume that the velocity, u, at St = 

0.5, has the form: 

ik.x. -ik-x^ 

u(x^,x^) = u(kj^,k 2 ) e e (2-28) 

If remeshing is done at St = 0.5, then, after remeshing, St = 
“0.5 and the new coordinate is 


x{ = Xj^ + X 2 


X.] = X2 


(2-29) 


Substituting (2-29) into (2-28) , we find the velocity after remeshing 
has the following form; 


^ ik.x -i(k^+k„)x 2 

u(xj|,xp = u(kj^ ,k 2 ) e e 


(2-30) 


If kji^ + k£ falls outside the available wavenumber range, the 
information will appear in lower wave numbers . This problem is called 
aliasing and can be eliminated by setting the amplitudes of wavenumbers 
at and above two-thirds of the maximum wavenumber to zero both before 
and after remeshing. 
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E . Computational Details 

In homogeneous shear flows, the length scales grow Indefinitely. 
Therefore, the assumption of periodic boundary conditions with finite 
period becomes Invalid after some time. However, If we use a relatively 
large computational box, we are able to study the flow for a limited 
time, l.e., until the large scales become comparable to the box length. 
The number of mesh points Is limited by the size of the computer memory. 
The maximum number of mesh points one can use on available computers 
Is 64 X 64 X 64 . In most of our computer runs , we were able to carry 
the computation up to St = 7.0, at which time the length scales become 
too large . The number of time steps for each run depends on the time 
step and the shear rate, and varies between 182 and 900. The computing 
time per computational time step Is about 68 .49 seconds of CPU time on 
ILLIAC IV. 

The acccuracy and stability of the numerical method are discussed 
In Appendix A and briefly explained below. A two-step method has two 
roots when applied to the equation; 

^ = (N^ + IN^) f (2-31) 
For the method of Eq . (2-11), these roots are 


1,2 


1 


/F 


tI (N + IN )^ - (N + IN ) (2-32) 

16 V c' ' V c“^ ^ ^ 


R 2 Is a parasitic root. For the case when f in Eq . (2-31) is Uj. , 
using linear stability analysis (see Appendix A), the Courant number, 
N^, , and the viscous stability parameter , N^ , are defined as 


N 

c 



Imax 


I + (1 + St) h 


2max 




At 


(2-33) 


N 


V 



^^2max 


'^Imax 

^11 

_ 


vAt (2-34) 
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The values of Rj^ and R 2 as a function of and are given In 
Table 2-1. We chose At such that and and therefore Rj^ and 
R 2 lie in the range where the solution is stable and more accurate. 
For the isotropic case, S = 0, we chose the time step. At = 0.21, so 
that both and are less than 0.3. For the homogeneous shear 
flow, the time step is also chosen to give a Courant number of less than 
0.3. 

The initial three-dimensional spectra of energy and scalar quantity 
are similar to those used for the isotropic homogeneous case, and are 
given in Fig. 2-4. 


20 



Table 2-1 


Accuracy and Stability Limits 


Nv 

Nc 

^1 

»2 

0.0 

0.0 

1.00 

0.00 

0.0 

0.1 

l . OO "^ 

0.050 

0.0 

0.2 

l . OO '*' 

0.010 

0.0 

0.3 

1.00+ 


0.0 

0.4 

1.01 

0.198 

0.1 

0.0 

0.905 

0.552 

0.1 

0.1 

0.904 

0.0782 

0.1 

0.2 

0.900 

0.124 

0.1 

0.3 

0.895 

0.177 

0.1 

0.4 

0.888 

0.232 

0.2 

0.0 

0.822 

0.122 

0.2 

0.1 

0.819 

0.136 

0.2 

0.2 

0.812 

0.174 

0.2 

0.3 

0.800 

0.225 

0.2 

0.4 

0.780 

0.287 

0.3 

0.0 

0.750 


0.3 

0.1 

0.747 

0.211 

0.3 

0.2 

0.738 

0.244 

0.3 

0.3 

0.721 

0.294 

0.3 

0.4 

0.694 

0.360 
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Fig. 2-1. The computational box. L|j is the box length and h is the mesh 
size. 


22 



WAVE NO. 


Fi«j. 2-2. Three-dimensional initial energy (or scalar) spectrum (32 riiesh 
points). 
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Fiy. 2-3. Coordinate transformation and rerneshing. 

(a) Coordinate at St = 0. 

(b) Coordinate at St = 0.5 (before remeshing). 

(c) Coordinate at St = -0.5 (after remeshing). 



icf 


1 1 1 1 t 1 1 — I — I 1 1 1 I I 

10 ' 

WAVE NO. 

Fig. k!-4. Three-diKiensional initial energy (or scalar) spectrum (64^ mesh 
points). 
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Chapter III 


COMPUTER PROGRAMMING AND CODE TESTING 


In this chapter we describe the construction and testing of the 
computer code and list the parameters of the flow fields simulated by 
the code.. 

3.1. ILLIAC-IV Computer 

The computations are done on the ILLIAC-IV computer. The ILLIAC 
contains 64 processors which operate in parallel and are controlled by a 
central managing computer. Each processor has a rather small memory. 
The array of processors has access to a large rotating disk on which 
most of the flow field data reside. Portions of the data are trans- 
ferred to the computer memory for the actual processing. The computer 
program and data management must be constructed carefully so that (a) 
the time spent transferring the data from the disk to the computer 
memory and vice-versa is minimized, (b) modification of the code and 
analyzing the results are easily done. We chose the data-management 
siystem developed by Lomax (1967). In that system, called the "pencil 
system," we move a column of data containing 8 x 16 x n words, where 
N is the number of mesh points in a given direction, at a time (see 
Fig. 3-1). All derivatives in one direction are computed while the 
appropriate column of data is in the computer memory. To compute all 
the derivatives, it is necessary to go through the entire flow field 
three times. This involves a lot of data transfer. However, the time 
consumed in transferring the data from the disk to the computer memory 
or vice-versa is a rather small portion of the total time (8.8% for 
64 mesh points). 

3.2. Tlie Code 

The computer program contains three different codes: 

A. A code which computes the initial conditions. This code gener- 
ates and stores the initial fields for both hydrodynamic and the scalar 
fields (see Chapter II, Section 2.3, for details of this program). 
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B. The main code, which advances the velocity and scalar fields in 
time using the Poisson, Navler-Stokes , and diffusion equations, as 
discussed in Chapter II. This code also stores the velocity and scalar 
fields on tape at the time when the coordinates are orthogonal (St = 
0); in the case of an isotropic homogeneous flow, the flow field is 
stored every 50 time steps. 

C. A code which uses the flow fields to compute various turbu- 
lence statistics and one- and three-dimensional spectra of various 
terms . 

We also used IBM 360 and CDC 7600 computers for final data reduc- 
tions and for plotting the data. 

The ILLIAC codes can operate with 16^, 32^, and 64^ mesh 

points. The run time for each case is given in Chapter II. 

3.3. Testing the Code 

It is important to make sure that the code is properly constructed. 
To check the validity of the results, the following items were examined. 

A. Conservation Laws 

The simulated flow fields satisfy continuity and the energy balance 
throughout the computations. The relative error for continuity is 

^'^i i^max 

— = 0 . 001 % 


and, for the energy balance, 

dq/2 

dt 

Production - Dissipation 


1 


few percent 


where the production and dissipation terms are calculated directly, and 

O 

a second-order approximation is used to calculate dq /dt. The errors 
in the energy balance are mainly due to the time differencing. 
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B. Comparison with Similar Computed Flow Fields . 

The simulated hydrodynamic field is checked against the results 
obtained by Rogallo (1978) and Feiereisen et al. (1981), who indepen- 
dently developed similar codes to simulate homogeneous turbulent shear 
flows for incompressible and compressible flows, respectively. The 
overall behavior of turbulent statistics for the hydrod 3 mamic field, as 
will be explained in Chapter VI, are similar in all three simulated flow 
fields . 


C . Comparison with the Exact Solutions of the Navier-Stokes and 
Diffusion Equations 


An exact analytical solution to the Navier-Stokes and diffusion 
equations for two-dimensional incompressible flow exists. This solu- 
tion, which is known as the vortex-cell or two-dimensional Taylor-Green 
solution, represents a flow field consisting of two-dimensional vorti- 
ces, arranged rectangularly, which decay in strength with time. The 
solution has the following form: 


-(kj+k^)vt 

u = - A cos k^Xj^ sin k^x^ e 

-(k^+k2)vt 

V = A sin kj^x^ cos k^-x.^ e 

^ 2 2 

✓ I 1 . ^ \ T<v._ 


-(kj^+k2)Dt 


C = B cos kj^Xj^ cos ^2^2 ® 

A^ *^1 -2(kj^+k2)vt 

- cos 2kj^Xj^ + — 2 " cos 2 k 2 X 2 e 

*‘2 


We set the initial conditions by putting t = 0 in Eqs (3-1) . We ob- 
tain the results for several values of k^^ and k 2 and for both 
inviscid and viscous flows. Note that k^^ and k 2 must be small to 
avoid aliasing. The exact results were compared with results obtained 
from the computer code after 80 time steps. To check the coordinate 
transformation, we obtained the results for skewed coordinates when St 
= 0.5. For all cases the results were obtained using 16 mesh points, 
and the difference between the exact results and the computer results 
was less than 1.7 x 10~^%. 
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D. Comparison with Experimental Data 


As mentioned in the first chapter, homogeneous shear flow was gen- 
erated experimentally by several people, we mainly check our results 
with the experiment by Tavoularis (1978) , because they generated both 
hydrodynamic and scalar fields. As will be discussed in Chapter VI, our 
results are in very good agreement with the results obtained from this 
experiment. However, since the experimental results are obtained at 
higher Reynolds number than the simulated results, there are slight 
differences between the results. 


3.4. Description and Details of the Computer Runs 

We have simulated several flow fields for both isotropic homoge- 
neous turbulence and homogeneous turbulent shear flow. For each run, 
the computed flow field can be divided into three parts. 


A. The "Developing" Region 


As described in the last chapter, we start the computations with 
artificial initial fields. Therefore, there is an initial period in 
which the simulated flow fields are not true turbulent fields. This is 
called the "developing" period. In the case of isotropic homogeneous 
flow, the development is considered complete when the velocity and the 
scalar derivative skewnesses reach asymptotic values. The skewness is 
defined as: 


S 


k 



(no summation) 


(3-2) 


In the shear flow with mixing of a scalar whose mean is linear in 
space, the criteria by which the end of the development period is judged 
are; 

(1) The shear stress correlation coefficient <Uj^U 2 >/u£u 2 and one 
point velocity scalar correlation coefficient <6uj^>/6'uj^ 
reach asymptotic values (see Fig. 6-17), where u| and 6' 
are the r.m.s. velocity and scalar fluctuations. 

(2) The three-dimensional spectra of kinetic energy and the scalar 
reach asymptotic shapes. 
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(3) The ratio of production to dissipation in the turbulent ki- 
netic energy and r.m.s. scalar fluctuation equations reach 
asymptotic values. 

(4) The rate of change of 8e5molds stress anisotropy tensor with 
time becomes small. 

B. The "Developed" Region 

After the initial period, there is a period in which the simulated 
flow fields can be regarded as true turbulent flows. In this period, 
the statistics are extracted from the flow fields. These will be 
discussed extensively in Chapters V and VI. For all computed runs, the 
developed region lies between St = 4 and St = 7 , where St is non- 
dimensional time. 

C . The "Anomalous" Regi on 

In both Isotropic and shear flows, the scales of motion grow with 
time, and there is a time at which the scales become too large for the 
computational box and the periodic boundary conditions are no longer 
accurate. Beyond this point, the flow fields do not accurately repre- 
sent homogeneous shear flow (or isotropic flow). Therefore the simu- 
lation is stopped. This time Is chosen by looking at the two-point 
correlations of the velocity and the scalar fields. The two-point 
velocity correlation is defined as 

^ < u. (x) u.(x+r) > 

(3-3) 

^ < U.(x) Uj(x) > 

3 

For the isotropic flows , we used 32 mesh points and various 

values of Reynolds and Prandtl numbers. The results will be discussed 

in Chapter V. 

3 

For homogeneous shear flow, we used 64 mesh points and simulated 
fourteen scalar fields and six different flow fields. Table 3-1 gives 
some details of the runs. 
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Table 3-1 


DETAILS FOR EACH SIMULATED FIELD 


Run 

* 

Reynolds 

Number 

Shear Ntunber 
SL/q 

Prandtl 

Number 

S3nnbol 

HS64A 

44.56 

18.55 

1.00 

□ 

HS64B 

19.00 

23.73 

1.00 

0 

HS64C 

19.00 

23.73 

5.00 

V 

HS64D 

44.56 

18.55 

5.00 


HS64E 

44.56 

18.55 

0.20 


HS64F 

19.00 

23.73 

0.20 

o 

HS64G 

60.35 

66.88 

1.00 

A 

HS64H 

60.35 

66.88 

5.00 

o 

HS64I 

60.35 

66.88 

0.20 

s 

HS64J 

24.20 

11.23 

1.00 

-f 

HS64M 

29.16 

8.94 

1.00 

X 

HS64N 

29.16 

8.94 

.20 

m 

HS640 

64.31 

6.38 

1.00 

O 

HS64P 

64.31 

6.38 

0.20 

® 


All values are given at St = 4.00. 
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Fig. 3-1, The computational box with three columns of 
data in different directions. 
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Chapter IV 

STATISTICAL MODELS OF TURBULENCE 


In this chapter, we present averaged equations for both the hy- 
drodynamic and scalar fields. We are interested in comparing and con- 
trasting the behavior of terms in the equations for two fields. We also 
present some of the models which have been used for the various terms in 
these equations. The terms which are modeled can be calculated from the 
numerical results of the preceding chapter and will be discussed in the 
following chapters. 

The equations are derived for homogeneous shear flow. They can be 
used for isotropic homogeneous flows by setting S = S2 “ In order 
to avoid complexity, we derive the equations for the case when the coor- 
dinates are orthogonal, i.e., St = 0. 


4 . 1 Equations for Turbulent Statistical Quantities 

The equations of motion for homogeneous turbulent shear are Eqs. 
(2-22) through (2-24). The momentum equation and the convection- 
diffusion equation for a scalar quantity have the following forms (note 
that we set St = 0): 


8u . 
1 

3t 


.1 - 


- p + VU . . 

i.JJ 


(4-1) 


and 


ae 

at 






(4-2) 


where, as before. 


S 



and 



Assuming the boundary conditions for two fields are similar, one sees 
that there is a close analogy between x^^— momentum uj^ and the scalar 
6 if ;Pr = v/D =1, S = S 2 » and ap/ax is small compared to the 
other terms. Similarly, if ap/3x2 or 3P/3x3 and S 2 U 2 are small 
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(i.e., in the isotropic case), Pr = 1 and 8=82, then an analogy 
exists between X2 or X3 -momentum and the scalar. 

The equations describing the three components of turbulent kinetic 
energy can be derived by taking the volume average of Eq. (4-2) multi- 
plied by Uj_ (no summation): 


8 < uj >/2 ^ 


A. 

at 


2 '*’11 

“ 2 ^11 ~ '^l'^2 ^ ^ 

(4-3) 

3 < »2 

>/2 

^ A 
2 ^22 

-is 

(4-4) 

at 


2 22 

»< “3 

>/2 

1 A 

2 ^33 

1 

— p 

(4-5) 

at 


2 33 


where 


♦ij ■ ^ 




(4-6) 

(4-7) 


and < > denotes the volume average. The third term on the right-hand 
side of Eq . (4-3) is the production term. Through the action of this 
term, energy is transferred from the mean flow to the fluctuating part 
of the flow. The pressure-velocity covariance terms, transfer 

energy from < u^ > to < u^ > and < U3 >. Thus we expect < 0, 

4>22 ^ ‘1’33 ^ Also, + <^22 ‘i’33 “ 

these terms neither produces nor destroys turbulent kinetic energy. 
Thus , these terms merely transfer energy from from one component to 
another. 8ince energy is transferred from the mean flow to the large 
scales of turbulence, we expect the production term and <t>^j to be most 
important at low wave numbers. is the dissipation term. Its 

trace, + £22 + ^33 > positive-definite, and is important at 

relatively high wave numbers. We shall show direct evidence of the 
spectra of these terms in the following chapter. 


The turbulent kinetic energy equation can be obtained by summing 
Eqs. (4-3), (4-4), and (4-5), and is 
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(4-8) 


3q^/2 

~5t~ 


- < 


> S - V < u. ,u. . > 

yj i >J 


The first term represents production and the second, dissipation. 

The equation for the r.m.s. scalar quantity can be derived from Eq. 
(4-2) and is: 


8 < >/2 

8 t 


- < “ 2 ® > ^2 


- D < 0 



(4-9) 


The first term on the right-hand side of Eqs. (4-9) is the production 
term, and the second is the dissipation term. 

The Reynolds stress < Uj^u 2 > and turbulent heat flux < 0Uj > 
equations can be derived in a manner similar to those used in deriving 
the above equations, and have the following forms: 


3t ^ 

-U 1 U 2 > = 

<u^>S- 

‘*’12 ^12 

(4-10) 

It < 

- 6 u 2 > = 

< 4 > ^2- 

^2 ^2 

(4-11) 

It < 

6 uj^ > = ■ 

- < U20 > S - 


(4-12) 


where defined by Eq . (4.6) and 


4>i = < P9,i > (^-13) 

and 

= (D+v) < > (4-14) 

In the above equations, the terms < u^ > S, <u| > S 2 , and - [< U 26 > S 
+ < Uj^U 2 > S 2 ] are production terms. 

It is not obvious whether £ 3^2 * ^2» ^1 h®lp prevent the 

growth of the correlations. However, from our results, they are nega- 
tive and are thus dissipative. The pressure-velocity and pressure- 
scalar covariances also tend to hinder the growth of the correlations. 
This can be shown by considering the return to isotropy at high Reynolds 
number. In this case the dissipation terms are negligible; thus the 
pressure-velocity and pressure-scalar covariance are responsible for 
destroying the correlations, and must therefore be negative. 
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The equations given above are exact, but they contain many more 

variables than there are equations, so the system is not closed. To 

close the system, it is necessary, to make assumptions about the rela- 

<3 OAAl-'i'rltS 

tionships among the various '-fiows. This is the issue of modeling that 
will be discussed next. 

4 .2 Statistical Models 

In this section we further consider the production, dissipation, 
and pressure-velocity and pressure-scalar covariance terms which were 
discussed in the last section. We review a number of models. 


A. Turbulent Reynolds Stress and Heat Flux 


The Reynolds stresses < > appear in the mean momentum 
equations , and the heat flux (mass flux if the scalar is concentration) 
< 6u^ > appears in the mean diffusion equation. If the only mean 
velocity gradient is 9Uj^/3x2, the mean turbulent Reynolds stress tensor 
takes the following form: 


< u,u. > = 

i 3 


< uj > 

< Uj^U^ > 


< Uj^U^ > 

< xxl > 

0 


0 

< u^ > 


(4-15) 


The principal stresses are 


2 2 

< u^ > + < U2 > 

- = + 


- ^^2 = 


2 2 
< u^ > + < u^ > 


2 2 2 

« u“ > - < U2 » 


2 2 2 
« u^ > - < U2 » ^ 


+ < UjUj > 


+ < > 


1/2 


1/2 


(4-16) 


and 


- Oo 


^3 > 
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The direction of the principal axis in the , X2 plane relative to 
the coordinate system is 


a 


1 



< uj > - < > 


(4-17) 


To close the mean momentum and diffusion equations , the Reynolds 
stress tensor and the heat flux vector must be modeled. The most common 
models for those terms are Newtonian-like constitutive equations, which 
is the Boussinesq eddy-viscosity concept. In the most general form the 
Reynolds stress tensor can be modeled by 


where 


< U^Uj > = 


^''lijpq ^pq 


1. 1 2. 

+ 3. q 6 


IJ 


(4-18) 


S 

pq 



u ) 

q»p 


(4-19) 


q 


2 


< Uj^Uj^ > 


(4-20) 


and '^T'ljpq ^ fourth-rank eddy viscosity tensor which needs to be 
defined. We looked at models for '^xijpq which are linear and quad- 
ratic functions of and ®'or each case the model was reduced 
to a simple expression which can be obtained using the following form of 
the eddy viscosity model: 


< . H- 1 (4-21) 

where is a scalar eddy viscosity. This model is not accurate for 

2 2 

this flow because it requires < Uj^ > = < U2 >♦ A more complicated 
model for < > was derived by Wilcox and Rubesin (1980). This 

model is especially good for the log region of the boundary layer when 
the magnitudes of production and dissipation are close. The model pre- 
dicts the behavior of < u^ > , < u^ > , and < u^ >. For Reynolds shear 
stress, < uj^\i2 >, the model reduces to Eq. (4-21). Rodi (1976) also 
derived an algebraic model for < u^Uj > which is nonlinear in < u^^Uj > 
and is based on the assumption that < u^uj >/< u^u^^ > = constant and is 
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good for high Reynolds number flows. The model (4-21) is used by many 
people to calculate the Reynolds shear stress < Uj^U 2 >• We shall use 
this model for < Uj^U 2 > and calculate the eddy viscosity. Modelers 
have used various forms for v^. We describe the ones applicable to our 
flow and specifically exclude models designed to apply only to boundary 
layers. 

on 

The simplest model for is based entirely pti information from 

the mean flow. 

v„ = C /S, .S . . (4-22) 

T ij xj m 

where C is a constant and is a mixing length. Models of this 

kind are used in boundary layers (for which is a function of the 

distance from the wall) and free shear layers (for which is a given 
fraction of the shear layer thickness). 

A more complex model for the eddy viscosity, which involves some 
information about the fluctuating field, has the following form: 

where C is constant and the turbulence velocity, q, is obtained by 
solving the turbulent kinetic energy equation, in which the higher-order 
terms are modeled (see W. C. Reynolds, 1976, for more details). is 

a mixing length which is prescribed as above. This is known as a one- 
equation model. 

A still more complex model is obtained by using Eq. (4-23) with a 
length scale defined by 

3 

L = (4-24) 

d o e 

in place of L^; here, e is the dissipation or, more correctly, the 
energy transfer from the large scales to the small ones , and can be ob- 
tained from an equation for the dissipation (see W. C. Reynolds, 1976). 
This is a two-equation model. For the low Reynolds numbers, the dissi- 
pation length scale is modified as follows: 
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3 C 

L , = C I- 1 + ^ 

d o e Re. 


(4-25) 


where Re is Reynolds number based on a Taylor microscale and q . 
Similarly, the heat flux vector can be modeled as follows; 


< 6u > = D -- — 

i ij 3x^ 


(4-26) 


where D^j is turbulent eddy diffusivity and is a second- rank tensor. 
If 8U^/8x^ is the only nonzero mean velocity gradient, has the 
following form; 


. = D 

ij 


ii 

^12 

21 

^22 


(4-27) 


D^j is not necessarily symmetric. To model this term, we decompose it 
as follows: 




(4-28) 


where 


“u ■ I- <“ij * “ji> - j 


(4-29) 


is a symmetric trace-free tensor , 


j (“ij - “ii> 


(4-30) 


is an antisymmetric tensor, and 


i D, , 6, . 

ij J kk ij 


(4-31) 


is a diagonal tensor. The above three expressions for need to be 

modeled,. In general, can be written as: 
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D. . = D, .(b. S. . , 5^. scalars) 

ij Ij^ Ij’ ij’ ij’ 


(4-32) 


where is the Reynolds stress anisotropy and is defined as 


< u .u . > , 

('•-w 


is defined by Eq. (4-19), and is the mean tensor defined as 

(4-34) 


i (U, . - U. .) 
ij 2 i,j j,i 


Note that we do not include the dissipation anisotropy, j » the 

right side of Eq. (4-32), because dj^j results mainly from the small 
scale structures, while D^j is dependent mainly on the large scales. 

and bji^j are both symmetric, trace-free, and have the same 


ij 


zero components, 
ity fluctuations 
velocity field. 

bi.. 


Since the scalar fluctuation is similar to the veloc- 


and 

ij 


b^. are both quadratic functions of the 

^ ( 1 ) 

Therefore, we assumed D;. is a linear function of 


= C b 
ij 1 S ij 


where Cj^ is a constant and L is a length scale 




is an anti- 

s)nnmetric tensor; therefore it must be modeled with the only antisymmet- 

d(2) 

ij 


ric tensor available, ^Ij* Neglecting higher-order terms, 
should have the following form: 

o 

.(2) . 


D 


ij 


2 


where C 2 is a constant. 

Finally, can be written as follows: 

n(3) = 


ij 




Therefore, Eq. (4-28) becomes 
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where 

numbers 


77s 




- + .C 2 g + •C3'Sij 

f 

W 13 


(4-33) 


2 ^, .C 2 , and .C 3 may be functions of^ Reynolds shear, Prandtl 
, and invariants of the above-mentioned tensors . 


B . Turbulent Prandtl Number 

The conventional turbulent Prandtl number (or turbulent Schmidt 
number if passive scalar is concentration) is defined as 

Pr„ = (4-36) 

^ r\ 


. . ac 




2 


/ft s 


The turbulent Prandtl number measures the degree of analogy between 
momentum and heat transport in turbulent flows and is required for heat 
transfer calculations . For flows involving mixing of chemical species , 
the important transport property is the diffusivity and the counterpart 
of the turbulent Prandtl number is turbulent Schmidt number. Reynolds 
(1885) proposed the Reynolds analogy, and Prandtl (1925) gave his 
mixing-length theory, both of which assume that Pr^ = 1.0. However, 
this assumption has been challenged by experimental results . Several 
models have been suggested (see, for example, A. J. Reynolds (1974)). 

We shall consider models which are recommended by several authors 
or used by many people in this field. The models we shall look at are: 

a) A. J. Reynolds proposed an empirical model based on experiments 
in wall flows, free turbulence, and core flows. Pr^ should depend on 
the molecular Prandtl number and a Reynolds number based on local turbu- 
lence intensity. The latter can be replaced by v.j/v, which is assumed 
to be a function of this Reynolds number. Based on these observations, 
he suggested that Pr,p be modeled as follows ; 
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Pr„ = C, exp - C,Pr 




(-) 


(4-37) 


where Cj^, C 2 , C^, and are positive constants. 


b) Wassel and Catton (1973) used mixing-length theory and assumed 
that eddies are spheres which move normal to the mean flow and have a 
radius proportional to the mixing length. Their model, based on this 


assumption, is 


-C^/CPr 

' %3/(v,7v) ' 
1 - e 


(4-38) 


Cj^, C 2 , and are parameters, and they suggest the values = 

1.05, C 2 = 5.25, and = 5. 

c) A model which simulates the idea that an eddy tries to equili- 
brate with the surrounding fluid has the following form: 


r 2 

a 

T 


+ aCj^Pe^ - (C^Pe^) (1-e 


-(a/C,Pe_) 


iPep)^ 


(4-39) 


where Pe^ is the turbulent Peclet number , 


and a and are constants. Crawford (1976) suggested a = 1.08 

and C = 0.2. 

d) Using the models for turbulent eddy dlffusivlty, Eq. (4-35), 
and turbulent eddy viscosity, Eq. (4-21), the turbulent Prandtl number 


becomes: 


< u: > 

S— -4 


(4-40) 


where and are constants. They may be functions of Reynolds 

number, shear number and Prandtl number, and Invariants of various ten- 
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C. Pressure-Strain Covariance 


The pressure-strain covariance, “l^i j » defined by (4-6). <()^j 

is responsible for exchanging energy between the various Reynolds 
stresses . 

The equation for the fluctuation pressure is obtained from the 
Navier-Stokes equation, as explained in Section 2-1. For the homogene- 
ous shear flow, the equation has the following form: 


P,li 


= - 2u 


4U. 




(4-41) 


The right-hand side of this equation contains two parts. The first part 
is proportional to the mean strain, j, and given a sudden change 

of Uj j, it changes rapidly. On the other hand, the second term on 

^ iJ 

the right of Eq. (4-41) does not depend explicitly on the mean strain, 
and it takes some time to change after a sudden change of the mean 
strain. Therefore, the two terms behave differently and need to be 


considered separately. We let p = 

p^^^ + p^^^ , where 


(1) 

P,ii 

- 2u U 

i.j j.i 

(4-42) 

p(2) = 

P,il 

- u. .u . . 
1.1 J.i 

(4-43) 


Following Lumley and Khajeh-Nouri (1974) , for homogeneous flows 
p^^^ can be defined in terms of the Fourier transform of the velocity 
field. By definition, 


p^^\x) 


/ 


p^^^(k) 


ikx j, 
e dk 


(4-44) 


In homogeneous flows, the mean velocity gradient is constant, so Eq . 
(4-42) gives 


2 

where k = . 

with p^^^ is 





(4-45) 


The part of the pressure-strain covariance associated 
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which can be reduced to (Lumley and Kha jeh-Nouri , 1974); 


= 2U, , ff -i u,(k) ru*(k')k’ + u*(k')k'l dk'dk (4-47) 

The double integral on the right-hand side of the above equation is a 
fourth-order tensor, so we can write 



2U G 

j,i ijpq 


(4-48) 


From continuity, ^Hpq °° also ®ijjq “ ^ “i“q ^^is tensor can- 
not be simplified and must be modeled. Following Lumley's argument, 
Gi jpq should be linear in the Reynolds stress , and it can also be a 
function of scalar invariants of the mean strain and Reynolds stress 
anisotropy tensor. The most general form of the pressure-strain term 
consistent with the above is : 


. I (1«, ) S ,2.3 
"^pq 5 1' pq^ 2 


R.S, +R,S, -i-<u.u.>U 

pk kq qk kp 3 i j 


where 


*i[' 

" T 1_T TT ^,k\pj 

S +U ^Is the mean strain, 

pq 2 \ p,q q,p/ 


R„^ = < u„u„ > is the Reynolds stress , 

pq p q ’ 


.6 

>J pq_ 


(4-49) 


(4-50) 


0 = V (0 ~ U ) is the mean rotation tensor, 

pq 2 ' p,q q,p" 


and is constant. The value suggested for this constant by dif- 

ferent workers varies between -1.23 and -2.5 (see W. C. Reynolds 
(1976) for more details) . 

The above model for our case, where j = 0» except for 

Uf 2 = S, has the following form: 



^11 

= < 


"1.1 

> 


TTAi ® < '^1^2 ^ 


'*’22 

= < 

pU) 

S,2 

> 

2 

3 

+ 1 A) S < > 

(4-51) 

S3 

= < 


'^3,3 

> 


S < U2.^2 ^ 



and 


♦12 ■ < 

+ - 15 - A^J S < > (4-52) 


( 2 ) 

The part of pressure-velocity covariance associated with p'^ ' is 




(4-53) 


and the total pressure-strain term is 


4, = .^(1) + <^(2) 


(4-54) 


Rotta (1951) proposed the following model for 


♦ij’ - -Vhj 


(4-55) 


where is the anlsotropis part of the Reynolds stress tensor and 

defined by Eq. (4-33), and ky is constant. 

Lumley (1978) suggested including the dissipation anisotropy in the 
model, so 





(4-56) 


Here 
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T 



Is the anisotropy of the dissipation tensor and 
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Is the dissipation. Using the dynamical equation for ^ , it can be 
shown that must be bigger than 2. The term acts to inter- 
change energy among components of turbulent kinetic energy, but neither 
creates nor destroys total energy. Lumley also included a second-order 
anisotropy tensor In the model; however, the model (4-56) seems to work 
quite well for homogeneous flows (see next chapter). Lumley (1978) also 
argues that the proportionality constant In the model (4-56) Is a 
function of Reynolds number and first-, second-and third-order invari- 
ants of the anisotropy tensor. 


D. Pressure-Scalar Covariance 

The pressure-scalar covariance, (j)^, is defined by Eq. (4-13). 
This term, which appears in the < 0u^ > equation, is responsible for 
destroying the velocity-scalar correlation. Following the procedure 
used in the last section, we divide into two terms. involves 

the mean shear explicitly, and is the remainder. 

Let us first consider From Eq. (4-13), we get 





> 


From Eq. (4-45), the Fourier transform of the pressure is 


(4-60) 


Let us multiply both sides 
(4-61) multiplied by 6 and 



k k 
q i 



(4-61) 


by 6 and add to the conjugate of Eq . 
integrate over Fourier space. 
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A A ^ 
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(4-62) 
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By definition, the integral on the left is the volume average of 
p .0. On the other hand, in a homogeneous flow, 

<p.0> = -<p6.> (4-63) 

,X ,1 


Using B;qs. 


(4-49) and (4-52), Eq. (4-51) becomes 


<!> 


( 1 ) _ 


2U 


’V J 


e dk 


(4-64) 


The integral on the right-hand side of (4-64) is a third-rank tensor. 
It cannot be computed analytically and needs to be modeled. Let us 
define as follows: 


where 



= - 2U 


q,j®qij 




dk 


(4-65) 


(4-66) 


A model for must satisfy the following conditions: 

• From continuity, Sqiq ~ 

• From the definition of g . g . = < 0u. >. 

®qij ^qqj j 

• From symmetry, g^j^j = g^q^. 

By analogy to the assumptions applied to the pressure-strain model, Eq. 
(4.49), we assume Sqij ® linear function of 0Uj^. The following 
model satisfies all of the above assumptions: 


2 

’qiJ 5 


6, < 0u. > - f 6, . < 0u > + 6 ^ < 0u 

iq j ij q qj 


Thus the model for the fast term, is 


^>] 


(4.67) 


= - ^ < U . > \ 8 , < 0u, > - i 6, . < 0u > + 6 . < 0u, > 1 

i 5 q , j |_ iq j 4 Ij q qj 1 J 


(4-68) 
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For our case, where ^ j = 0, except for ^2 ~ S, the model has the 
following form : j 

< > = 2SAq < 0U2 > (4-69) 

. SA 

< -^6 2 ^ > (4-70) 

where A^ = 2/5. Comparing the above equations with Eq . (4-40), we see 
great similarity between the models for the hydrodynamic field and the 
corresponding models for the scalar field. 


Lumley (1978) argues that, as < 6uj^ > •> 0 , 0. He construc- 

ted an expression for 4’^^^ which satisfies all the above conditions: 


= - 2A U , 

1 o q,j 


qi 


< 6u > < 0u. > 

q L_ 

< 0u > < 0u > 


< 0u. > 
J 


(4-71) 


where A^ = 1/2. 

For our case , this model becomes 

< p (')0 


1 > = 2A S 

,1 o 


< 6u2 


< 0u^ > + < 0U2 > 


^ < 0ui > 


(4-72) 


< 0u, > < 0u. > 

< P 0 o > = - 2A S ^ = 5 - < 0u„ > (4-73) 

’ ° < 0u^ > + < 0U2 > 

Now let us consider the second part of the pressure-scalar covari- 


ance , 


(j,p) = < p(2)e^^ > 


(4-74) 


This term appears in the heat flux equation and tends to destroy the 
velocity-scalar correlation. The average heat flux equation in the 
absence of a mean scalar gradient is : 


3 < u^0 > 
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- < 


.(2), 


> - (v+D) < 0 u. > 

>J >J 


(4-75) 
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Both terms on the right side of the above equation tend to destroy the 
< Uj6 > correlation. To model these terms, let us define by 


- < 6uj > ^ = < > - < (vbD) 


Lumley argues that the principal axes of are aligned with the 

principal axes of j , where 


< U^Uj > < 0U^ > < 0Uj > 


ij 


(4-77) 


2 

0 q 


Therefore, must be a function of 

Assuming a linear relation, we get 


♦ij ■ Vij + Vij 


(4-78) 


where Aj^ and A 2 are functions of anisotropy invariance. 

Assuming A 2 = 0, the model for our flow will have the following 

form: 


and 


< P^^^0 -^ > - < 0 j^u^ j^> == .^< 0Uj^ > (4-79) 

’ . > J- 


< 2 > - < (^D) ® k“2 k ^ " ^1 ^ 

’ > » 


The above model in fact is very similar to Lumley 's model for the Rotta 
term, Eq. (4-56). 
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Chapter V 


DECAY OF A PASSIVE SCALAR IN HOMOGENEOUS ISOTROPIC TURBULENT FLOWS 


In this chapter we present results for decay of the velocity field 
and a passive scalar in Isotropic homogeneous flows. The decay rates, 

various length scales, two-point correlations, and one- and three- 
dimensional spectra of various terms are presented and discussed for 
both the hydrodynamic and scalar fields. The effects of the Reynolds 
and Prandtl numbers on the decay rates, microscales, and skewness are 
also considered. All of the results presented in this chapter were ob- 
tained using 32 X 32 X 32 mesh points , so the Reynolds number range is 
limited to Re^ < 25. 

5 . 1 Results from a Typical Run 

In this section, we present the results obtained for a typical run 
(IH32-2). This run has the following initial values: 

qX 

Reynolds number Re^^ = ~ = 11.00. 

Integral length scale ratio L/Lq = 1.30. 

Taylor micro-scale ratio X/X0 = 1.23. 

Integral length scale L = 0.82. 

Taylor micro-scale X = 0.23. 

Time step At = 0.0210. 

Computational box length L^ = 6.28. 

Initial turbulence velocity q = 0.02. 

Molecular Prandtl/Schmidt number Pr = v/D = 1.00. 

The results are shown in Figs. 5-1 through 5-13 and discussed 
below. 

A . The Ilydrodynamic Field 

Figure 5-1 shows the evolution of the three-dimensional energy 
spectrtam. The 3-D energy spectrum is obtained by summing the energy in 
spherical shells of inner radius k and outer radius k + Ak , where 


50 



k is the magnitude of the wavevector and Ak is the difference between 
the magnitudes of the two nearest wavenumbers. As stated in Chapter II, 
the Initial 3-D spectrum Is a square wave. The 3-D spectrum evolves 
from the artificial square wave Inertial condition to a realistic low 
Rejmolds number energy spectrum. Since the Re 3 molds number Is small, 
there is no Inertial subrange. At low wavenumbers, the slope becomes 
4 at long times — an expected result. The energy at wave numbers higher 
than 2/3 of the maximum wavenumber Is set to zero to avoid aliasing. 

Figure 5-2 shows the 3-D spectra of the components of the turbulent 
kinetic energy at the final time step. As can be seen, the flow is 
slightly anisotropic at low wavenumbers and Isotropic at high wavenum- 
bers. This, too, is expected, especially as there are only a few low 
wavenumbers . 

The lateral and longitudinal one-dimensional spectra of the veloc- 
ity field are shown In Fig. 5-3. The 1-D spectrum of the velocity field 
In the kj^ direction Is defined as 

EijCkf) = u*(k) dk 2 dk 3 (5-1) 

As shown, velocity field is isotropic, £ 22 (^ 2 ) “ expected. 

Figure 5-4 shows the time evolution of the turbulent kinetic energy 
and Its three components In log-log coordinates. As can be seen, ne- 
glecting the "developing" region (defined In Chapter II), the turbulent 
kinetic energy is a straight line. If e(t) is the turbulent kinetic 
energy at time t , then 

e(t) = A t~“ (5-2) 

where n is the decay exponent. In this particular run, n = 2.4. The 
decay of the turbulent kinetic energy has been observed In many experi- 
ments. It has been shown analytically that the exponent, n, for very 
low Reynolds numbers Is 2.5. It has been shown both experimentally and 
analytically that the decay exponent for high Reynolds numbers is 1.20. 
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Figures 5-5 show the lateral and longitudinal two-point correla- 
tions at initial time and the final time. The. two-point correlations of 
the velocity are defined as 


R^j(r) = < u^(^r) > / < u^(x) u^(x) > (5-3) 


As shown in Fig. 5-5a, the initial two-point correlations oscillate 
sinusoidally. This is due to the narrow band nature of the initial 
spectrum. This effect vanishes with time. 


The integral length scale, L^jCr^), is defined as twice the dis- 
tance at which first reaches 0.1. This differs from the 
standard definition of the integral length scale, which is the integral 
of the two-point correlations. However, since the two-point correlation 
can have negative values at large r, the standard integral length 
scale may behave poorly. The time evolution of the integral length 
scales is shown in Fig. 5-6. 


The Taylor micro-scales are shown in Fig. 5-7. The Taylor micro- 
scale, is defined as the inverse of curvature of a two-point 
correlation at r = 0. Both the integral and Taylor micro-scales in- 
crease with time. However, the integral length scales increase faster. 


B. The Scalar Field 

Figure 5-8 shows the time evolution of the 3-D spectrum of the sca- 
lar quantity. It is very similar to the 3-D energy spectrum. However, 
the scalar fluctuation spectrum is slightly higher at high wavenumbers. 
This is also observed experimentally for the flows with Prandtl number 
close to 1. It suggests that the dissipation in the velocity field is 
lower than dissipation in the scalar field. 

The one-dimensional spectrum, two-point correlations, integral- and 
micro-scales of the scalar quantity are also similar to those for the 
velocity field. They are shown in Figs. 5-4 through 5-7, along with 
corresponding values for the velocity fields. 


52 



The decay of the scalar intensity < 0^ >, with time is shown in 
Fig. 5-9. As can be seen, its behavior is similar to that of the tur- 
bulent kinetic energy history and can be fit by: 

< 6 ^ > = B t"“ (5-4) 

where m is the decay exponent of the scalar and B is a constant. In 
this particular run, m = 2.25, which is within the range obtained from 
the experimental results. The behavior of <0 >, Eq. (5-4), is also 

observed experimentally. The magnitude of m obtained from several ex- 
perimental results (see Table 1-1, also Lumley and Warhaft (1978), and 
Domis (1979) for more details) varies between 0.87 to 3.0, depending on 
the initial length/scale ratio, l/l 0 , and scalar intensity; the Prandtl 
number was 0.7. There are no experimental results for m at low Rey- 
nolds numbers or other Prandtl numbers. 

^•2 Description of the Simulations 

Nine simulations of isotropic homogeneous flow were made. The 
simulations were obtained for various Prandtl numbers and a limited 
range of Reynolds numbers. The initial micro- and integral-length/scale 
ratios for all of the runs are fixed: = 1.23 and L/Lq = 1.30. 

Table 5-1 shows some statistics for the various runs. In this table, 
Re^ in the third column is the Reynolds number at the final time step. 

R is the ratio of the decay exponents: 

R = m/n (5-5) 

The decay exponent ratio, R, is an important parameter. It is 
the ratio of characteristic time scales of the velocity and scalar 
fields: 

R = — (5-6) 

< ® >/^e 

where e is the rate of dissipation of kinetic energy and £9 is half 
the rate of the dissipation of the scalar intensity, and is defined as 
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(5-7) 


= fD< 


0 .6 . > 


In the following section we shall discuss the effects of Reynolds 
and Prandtl numbers on the decay-rate exponents. 


5 . 3 Effects of Reynolds and P randtl Numbers on Some Turbulence Statis- 
tics 

As explained in the previous section, flow fields at various Rey- 
nolds and Prandtl numbers were generated. In this section we examine 
the behavior of the various turbulence statistics. 


A. The Length Scale Ratio 

Figure 5-10 shows the behavior of the microscale ratio, 1/^9, as 
a function of Reynolds nximber for various Prandtl numbers. As can be 
seen, the microscale ratio increases with Rej^ when Pr < 1, and it 
decreases with Re^^ when Pr > 1. On the other hand, at high Reynolds 
numbers, the length scale ratio is independent of Prandtl numbers. 
Ironically, as Rej^ •*■0, the length scale ratio reaches a constant. 
This constant depends on the Prandtl number. 

B. The Decay Rate 

Figure 5-11 shows the decay exponent n vs. Reynolds number. As 
can be seen, at very low Reynolds numbers n reaches 2.5. This has 
been shown analytically (see Batchelor (1948)). As the Reynolds number 
increases , the decay exponent n decreases . At high Reynolds numbers , 
however, it reaches an asymptotic value. W. C. Reynolds (1976) argues 
that the asymptotic value is 1.2. The experimental values for n at 
high Reynolds numbers, obtained by Warhaft and Lumley (1978) and Comte- 
Bellot and Corrsln (1971), is 1.25 ± 0.06. As shown in Fig. 5-11, the 
experimental value for n obtained from Tavoularis et al. (1978) is 
very close to our simulated results. The Bennett and Corrsin (1978) 
result, however, does not agree with our results. This may be due to 
the difference in origin of the decay rate. 
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C . The De6ay Rate Exponent Ratio 

The decay exponent ratio, R, defined by Eq. (5-5), Is shown In 
Fig. 5-12 as a function of Reynolds number for various Prandtl numbers. 
As can be seen from the figure, R increases with Reynolds number at 
Pr < 1, and it decreases with Reynolds number at Pr > 1. At high 
Reynolds number, R asymptotically reaches a constant. This constant 
depends on the Prandtl number. As Re^ 0, however, it seems that the 
decay exponent ratio reaches a constant. This constant is Independent 
of Prandtl number. 

D. The Velocity Derivative of Skewness 

Velocity derivative skewness is defined by Eq. (3-2). Figure 5-13 
shows the skewness vs. Reynolds number. The results shown in this fig- 
ure were obtained from a 16 x 16 x 16 mesh point code and various 
numerical and experimental results. As can be seen, skewness has a 
maximum value of approximately 0,5 at Re^ = 20. For Re^ > 20, 
skewness decreases and it seems to reach a constant at about Re^ = 
100. At low Re 3 molds numbers (l^^X ^ 20), the skewness decreases and 
Sk -»• 0 as Re^^ 0. 
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Table 5-1 


DESCRIPTION OF COMPUTER RUNS FOR DECAY OF A PASSIVE SCALAR 
IN ISOTROPIC HOMOGENEOUS TURBULENT FLOWS 


Simulations 

<^"a)o 


Pr 

^u/^e* 

n 

m 

r 

IH32--2A 

11.00 

3.00 

0.20 

0.634 

2.40 

3.00 

1.25 

IH32-3A 

22.30 

10.00 

0.20 

0.703 

1.70 

3.00 

1.76 

IH32-4A 

44.40 

22.00 

0.20 

0.830 

1.50 

2.90 

1.93 

IH32-2 

11.00 

3.00 

1.00 

1.27 

2.40 

2.25 

0.94 

IH32-3 

22.30 

10.00 

1.00 

1.23 

1.70 

1.70 

1.00 

IH32-4 

44.40 

22.00 

1.00 

1.19 

1.50 

1.61 

1.07 

IH42-2B 

11.00 

3.00 

5.00 

1.68 

2.40 

0.75 

0.31 

IH32-3B 

22.30 

10.00 

5.00 

1.40 

1.70 

0.50 

0.29 

IH32-4B 

44.40 

22.00 

5.00 

1.26 

1.50 

0.34 

0.23 


(Re^)jj is the initial turbulent Reynolds number. 

The initial microscale ratio for all of the flow fields is the same , 
and its value is (^u^^ 6 ^o ~ 1*23. The initial lateral integral length 
scale ratio for all flow fields is (Lu/L 0 )^ = 1.30. 
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WAVE NO, 

Fig. &-2. Three-dinenslona) spectra of the turbulent kinetic energy compo- 
nents at St ■■ 4.0. 
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Fig. S-3. One-dinenslotMl Sfiectr* of the velocity components end scalar 
fluctuations at St * 4.0. 
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Fig. S.4. Decay of the turbulent kinetic energy and its components with time. 
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Fiu. Taylor microscale ratio as a function of Reynolds number for 

various Prandtl niiiibers* 
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Fiy. 5-11. Decay exponent of turbulent kinetic eneryy as a function of 
Reynolds number. 



Fig. 5-12. Decay exponent ratio as a function of Reynolds number for various 
Prandtl nuiiibers. 
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Fig. 5-13. Velocity skewness as a function of Reynolds number 


Chapter VI 


RESULTS FOR HOMOGENEOUS TURBULENT SHEAR FLOWS 
WITH A MEAN GRADIENT OF A PASSIVE SCALAR 

In this chapter the results for both the velocity and the scalar 
flows in a homogeneous turbulent shear flow with a passive scalar, whose 
mean is linear, are presented and discussed. Six velocity and fourteen 
scalar fields were generated. The flow fields are obtained for various 
Reynolds, Prandtl, and shear numbers. In this chapter a complete set of 
results for one of the simulated flows is presented. Also, the effects 
of the parameters on various turbulence statistics are discussed and 
models for Reynolds stress, heat/mass flux, pressure-strain, and 
pressure-scalar covariances are tested and their validity is discussed. 

6.1 Dimensional and Non-Dimensional Parameters 

The important dimensional parameters of this flow are presented in 

Tables 6-1 and 6-2. The dimension of the scalar, 0, in Table 6-2, is 

(a) degrees if the scalar is temperature and (b) mass if the scalar is 

species. The integral and microscales are defined in the previous 

chapter. 

From the parameters in the tables , the following nondimensional 
parameters can be formed; 

Re = q^/v, Reynolds number. 

^< 4 -^ = SL/q, shear muaber. 

Pr = v/D, Prandtl number. 

= 6 '/S 2 L, nondimensional scalar fluctuations. 

L/L 0 = length scale ratio. 

St =■ nondimensional time. 

Note that L/A depends on Reynolds number and therefore is not 
included above. The shear number can be written as 

^ ( 6 - 1 ) 
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which is the ratio of the time scales of the large scale structures and 
the mean flow. Assuming L = q^/e, can be written in the following 
form: 





(6-2) 


which is proportional to the ratio of the production to dissipation of 
turbulent kinetic energy. 

The parameter is analogous to the Inverse of It is 
approximately the ratio of the production to dissipation of the scalar 
fluctuations . 


Some of the nondlmensional parameters mentioned above are related 
through the dynamics of the flow, that is, they are not independent. 
The results obtained from a set of 32 x 32 x 32 mesh point computer 
runs show that the length scale ratio, L/L 0 , and the nondlmensional 
scalar quantity, '^ 2 * dependent parameters. Figure 6-1 shows 

^ function of time for four different runs. The initial Rey- 
nolds, shear, and Prandtl numbers and length scale ratios are the same 
for all runs, but the initial is changed from one run to another. 

As shown in the figure, the final value of independent of its 

initial value. In other words, ''^^2 ® function of Reynolds, shear, 

and Prandtl numbers, i.e., a dependent parameter. 

Similarly, as shown in Fig. 6-2, the final value of the length 
scale ratio L/Lq is independent of its initial value, i.e., L/L 0 is 
a dependent parameter. Therefore, the only nondlmensional parameters we 
need to deal with are Reynolds, shear, and Prandtl numbers. The scalar 
field is determined by the velocity field and the Prandtl number. 


6.2 Description of the Simulations 

Six simulations of the hydrodynamic field and fourteen simulations 
of the scalar field have been performed. All of the simulations used 
64 X 64 X 64 mesh points. The results presented in the rest of this 
chapter are from these simulations. 
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As discussed in Chapter III, the simulations use an artificial ini- 
tial field. There is also a time at which the turbulent structures 
become large compared to the size of the computational domain, and the 
periodic boundary conditions no longer hold. The results are not accep- 
table after that time. According to the criteria mentioned in Chapter 
III, the period in which the data are acceptable lies between St = 4 
and St = 7 for all of the simulations. 

The simulated flow fields are stored on tape, and the turbulence 
statistics are calculated only when St is an integer. This is done 
mainly because, when St is an Integer, the coordinate is Cartesian and 
the calculations of the turbulence statistics involving derivatives are 
easier. 

Thus, four flow fields (at St = 4, 5, 6, and 7) are available 
from each simulation. Consequently, a total of 24 hydrodynamic fields 
and 56 scalar fields are available. The available flow fields cover a 
rather large range of Reynolds, shear, and Prandtl numbers. Table 6-3 
shows the range of these nondimensional parameters in our simulations 
and from the available experimental results. 

As can be seen from the table, the range of Reynolds number covered 
is lower than the available experimental values. We are limited to low 
Rej^ because there is a maximum Re^ beyond which we are not able to 
resolve the full range of scales of turbulence. As will be shown in the 
following sections, we are mainly interested in turbulence statistics 
which depend on the large structures and which may be only weakly Rey- 
nolds number dependent, so the Reynolds number limitation may not be 
very serious. This will become clear in later sections, where we com- 
pare our low Reynolds number results with the experimental results. 

As can be seen from Table 6-3, we have covered a rather large range 
of the shear number. Near the upper limit the nonlinear and viscous 
terms in the Navier-Stokes equations (Eq. (2-18)) become negligible, and 
one can apply rapid distortion theory. At lower values of the shear 
number, the flow becomes isotropic (the results for this case are dis- 
cussed in the previous chapter) . 
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The only available experimental data for the scalar field , obtained 
for Pr = 0.7, were taken by Tavoularis and Corrsin (1978). We have 
covered a significant range of Prandtl number (from 0.2 to 5.0). 


At very high Prandtl numbers (Pr » 1) , the diffusion equation 
(2-20) becomes; 


8t 


(6-4) 


Therefore, the results are independent of Prandtl number (we assume 
viscosity is fixed). On the other hand, for very low Prandtl numbers 
(Pr « 1) , the effect of shear becomes negligible and the diffusion 
equation is 


86 _ 8^6 

8t 8x.8x. 

1 J 


(6-5) 


and the results are discussed in the previous chapter. 

For the region between these two extremes , where Pr ~ 1 , there is 
no analytical solution. The effect of the Prandtl number on the various 
turbulence statistics in this region is studied and discussed , using our 
simulated results, in the following sections. 

Table 6-4 shows the nondimensional parameters for each simulation 
in the "developed" region. In Appendix C, we have provided a set of 
tables which give the raw data and various turbulence statistics from 
each simulation. These tables may be helpful to the readers in further 
investigations . 


6.3 Complete Results for a Typical Simulation 

In this section, we present the results obtained for a typical run 
(Run HS64J) , which has the following initial values: 

Reynolds number Re^ = 20.00. 

Shear number 3.30. 

Nondimensional scalar fluctuation x«^02 = 3.20. 

Prandtl number Pr = 1.00. 
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Integral length scale ratio L/L 0 = 1.40. 

Taylor microscale ratio = 1.29. 

Integral length scale L = 0.71. 

Taylor micro-scale X = 0.22. 

The above values are non-dimensionalized , as explained in Section 
2.3. The results will be discussed for both hydrodynamic and scalar 
fields in the following sections. 

A. The Hydrodynamic Field 

A-1. Three-Dimensional Energy Spectrum 

Figure 6-3 shows the evolution of the three-dimensional energy 
spectrum. The initial 3-D energy spectrum is a square wave. At a later 
time, the 3-D energy spectrum has a slope of 4 at low waventimbers , as 
expected. No inertial subrange with slope of -5/3 appears in the 3-D 
energy spectrum because the Reynolds number is low. Figure 6-4 shows 
the 3-D spectra of the three components of the turbulent kinetic energy 

at St = 4. As can be seen, there are anisotropy effects at low and 

2 

moderate wavenumbers. The 3-D spectrum of u„ has more energy at low 

^ 2 

wavenumbers and less at moderate wavenumbers , whereas u^^ has more 
energy at moderate wavenumbers and less at low wavenumbers. These are 
effects of the applied shear. 

A-2. Turbulent Kinetic Energy History 

The time evolution of turbulent kinetic energy and its three compo- 
nents are shown in Fig. 6-5. At St = 0, the flow field is nearly iso- 
tropic. At later times, however, the flow becomes anisotropic. The 
streamwise component of the kinetic energy becomes the largest while the 
normal component is the smallest. This is also observed experimentally 
in both homogeneous shear flows (Champagne et al. (1970) and Harris et 
al. (1976)) and boundary layers (Klebanoff (1955), Bradshaw (1967)). 

As shown in Fig. 6-5, the turbulent kinetic energy first decreases 
with time and then increases. This is because the production of turbu- 
lent kinetic energy is small at early times. The production increases 
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with time, and there is a time after which the production is larger than 
the dissipation and the energy begins to grow. This effect can be seen 
from the experimental results obtained by Harris et al. (1976) and 
Tavoularis (1979). The experimental results obtained by Champagne et 
al. (1970) and Rose (1970), however, do not show the increase in the 
turbulent kinetic energy. Their measurements did not go far enough 
downstream. 

From our computed results and the experimental results of Tavou- 
laris (1979) , the flow fields seem to be going toward some kind of 
structural equilibrium. Statistical equilibrium is a state in which the 
components of the Reynolds stress anisotropy tensor, reach 

asymptotic values. As noted before, we were able to carry out the simu- 
lations only up to St = 7 and did not reach an equilibrium region. 
Rogallo (1981) simulated this flow with 128 x 128 x 128 mesh points. 
He was able to continue the computation to St = 16 and found that the 
rate of change of components of b^^j with time decreases , but it does 
not become zero. 


Length Scales 


Figures 6-6 and 6-7 show the evolution of integral and Taylor mi- 
croscales in the streamwlse direction, respectively. The length scales 
were obtained from the two-point correlations in the manner explained 
in the previous chapter. Figure 6-8 shows the Kolmogorov scale, the 
mixing-length scales, and the transfer (or dissipation) length scale, 
where the Kolmogorov length scale, n, is defined by Eq. (2-1) and the 
mixing length is defined as 



and the transfer length scale, is defined as 


( 6 - 6 ) 



(6-7) 


As can be seen, all length scales Increase with time. As shown in the 
figure, the mixing length is nearly constant, while and n 
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increase with time. The Kolmogorov length scale is the smallest, and 
transfer length scale is the largest length scale. Thus the mixing 
length cannot be correlated with any of the other length scales. 

A"4 . Two-Point Correlations 

Figures 6-9 through 6-11 show the two-point correlations of the 
velocity (Eq. (5-3)) in the streamwlse direction at three different 
times. As can be seen, the two-point velocity correlations, R^j(r), 
at St = 0, are nearly zero at large r, i.e., there exist only small 
eddies in the computational box. However, at later times, the results 
show that the velocities become well correlated in the computational 
box. As discussed before, there is a time after which the large eddies 

become comparable with the box length (in other words, R^j(r) becomes 
large at large r) , and the computation has to be stopped. In most of 

our results, R. .(r) is less than 0.3 at r = L./2 at the end of the 
ij b 

computation, where is the computational box length. 

A-5 . One-Dimensional Spectra 

Figures 6-12 through 6-14 show the 1-D spectra of the velocity 
components in three directions. The difference between the three com- 
ponents of the velocity spectra show the anisotropy at various eddy 
sizes. As can be seen, even at high wavenumbers the three components 
are different, i.e., the small scales are anisotropic. This may be 
because the Reynolds number is low. The anisotropy at high wavenumbers 
will be discussed in more detail in Section 6-5D. The fluctuations of 
the spectra at the low wavenumbers are due to the small number of sam- 
ples at which spectra are calculated. 

A- 6 . Three-Dimensional Pressure Spectrum 

As discussed in Section 4-2C , the pressure consists of two parts, 
p^^^ is proportional to the mean shear, while p^^^ does not depend 
explicitly on the mean flow. Figures 6-15 and 6-16 show the 3-D spectra 
of each component, as well as the 3-D spectrum of the total pressure at 
various times. Understanding the behavior of each component helps in 
modeling; the pressure-strain and pressure-scalar covariances. As shown 


71 


In Fig. 6-16, is larger at low wavenumbers, i.e., it depends 

mainly on the large scale structures. On the other hand, p^^^ is 
larger at high wavenumbers, i.e., it depends mainly on the small scale 
structures. We also see that the 3-D spectrum of p^^^ has a slope 
of -5 at high wavenumbers. This is observed in all of the simulated 
flow fields. Also, the k~^ part of the spectrum is observed to become 
more important with time; the peak in the spectrum moves upward and to 
the left as time goes on. 

A-7 . One-Point Correlations 

The one-point correlations of the velocities are defined as 

< u^u . > 

C. . = — I — — (no summation) (6-8) 

ij 1 j 

Note that = 1. Figure 6-17 shows the one-point velocity correla- 

tions as a function of time. The only non-zero component for (1 ^ j) 
is < >/u’u 2 , which starts from zero and reaches an asymptotic 

value of about -0.6. Similar values were obtained from computed flow 
fields by Feiereisen (1981) and Rogallo (1981). The experimental value 
lies between -0.4 and -0.5. The reason for this difference may be 
the Reynolds number difference between the computed and the experimental 
results . 

The one-point correlation, Cj ^2 Reynolds number dependent, be- 
cause < Uj^U 2 > depends mainly on the large scale structures and: is 
therefore not strongly Reynolds number dependent, whereas u| and 
U 2 depend on both the large and small scales and are thus more Rey- 
nolds number dependent. Therefore, as the Reynolds number decreases, 

< Uj^U 2 >/uJ^U 2 increases. The effect of Reynolds number on this quan- 
tity is given later in this chapter. 

A-8. Reynolds Stress Anisotropy Tensor 

The Reynolds stress anisotropy tensor, ^i j > defined by Eq. 

(4-44). Figure 6-18 shows b^^j as a function of time. If the compo- 
nents of b^j reach asymptotic values, the flow has reached structural 
equilibrium. As shown in the figure, bj^j start from zero, indicating 
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that the initial velocity field is isotropic, and change with time. At 
early times the change rapidly; later, however, the variation of 

is smaller. It seems that the flow might reach a state of equilib- 
rium if the computation were continued. 


A-9. Dynamical Reynolds Stress Equations 

The volume-averaged equations for the Reynolds stress , < ' Uj_Uj > , 

are given by Eqs. (4-3)- (4-5) and (4-10) and discussed in Chapter IV. 
The equations contain production, pressure-strain, and dissipation 
terms. Transfer terms do not appear in these equations because the vol- 
ume averages of these terms are zero. Figures 6-19 through 6-22 show 
the values of each term in the various stress equations as a function of 

time. As can be seen, the pressure-strain term is negative in the 
2 

< u^ > equation and positive in the equations for the other components 

^ 2 
of the Reynolds stress. That is, energy is transferred from < u^^ > to 

the other components of < n^uj > through the action of the pressure- 

strain terms. 


Figure 6-23 shows the terms in the volume-averaged turbulent ki- 
netic energy equation (4-8). No pressure-strain term appears in this 

3Ui 

figure because, by continuity, < p > is zero. 

i 

Figures 6-24 through 6-28 show three-dimensional spectra of various 
terms in < > equations at St = 4. The production and pressure- 

strain terms are concentrated on the large scales , whereas the dissipa- 
tion terms mainly come from the small scales, as expected. As can be 
seen in Fig. 6-25, unlike its behavior in the other equations, the 

^“2 2 

pressure-strain term, < p ^ — > , in the < u„ > equation is not only 

^^2 ^ 8 u 2 

a source term but it acts as a transfer term, too. That is, < p > 

2 

takes energy from moderate-sized eddies and transfers it to the small 
and large eddies. 


The transfer terms have the following form: 


T 


ij 




+ u, 
J 




(6-9) 


73 



The 3-D spectrum of T. . is defined by 

A J A A 



( 6 - 10 ) 


As can be seen from Figs. 6-24 through 6-28, the transfer term transfers 

energy from both small and large scales to the medivun-sized eddies, 

2 

except for the transfer term in the < u^ > equation, which transfers 
energy from large eddies to the small eddies. These effects, which are 
somewhat surprising, have been observed in all the calculated flow 
fields. McMillan and Ferzlger (1978) have also found that the transfer 
term is negative at high wavenumbers when the shear is high. These 
results were obtained at high Reynolds numbers. 


A-10. The Dissipation Anisotropy Tensor 

The dissipation anisotropy tensor, ‘if j > is defined by Eq. (4-47). 
Figure 6-29 shows d^j as a function of time. It has generally been 
assumed that the small scales in turbulent flows are isotropic, i.e., 
the dissipation anisotropy is nearly zero. However, as can be seen, the 
components of d^j start from zero and vary with time, dj^j^ increases 
while d 22 decreases in the developing region. It also can be seen 
that d ^3 is an order of magnitude smaller than the other components 
and is therefore not significant. Furthermore, as we shall see later, 
d^j seems not to be very Reynolds number dependent, and this result is 
of considerable significance in turbulence modeling. Also note that the 
dissipation anisotropy is very similar to the Reynold stress anisotropy 
(Fig. 6.18); its components behave in the same way with time. 


B. The Scalar Field 

The length scales, two point-correlations, and one-dimensional 
spectra of the scalar fluctuations behave in a manner similar to those 
of the velocity field. They are shown in Figs. 6-7 through 6-14. From 
the two-point correlations and the 1-D spectra of the scalar and veloc- 
ity fields, it can be seen that the scalar quantity behaves very much 
like the streamwise velocity. 
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The time evolution of the 3-D spectrum of the scalar fluctuations 
is shown in Fig. 6-30. As can be seen from the figure, the behavior of 
the 3-D spectrum of the scalar quantity is similar to that of the 3-D 
energy spectrum. The scalar-fluctuation intensity as a function of time 
is shovm in Fig. 6-31. The scalar intensity Increases with time in the 

developed region. Its behavior with time is similar to the behavior of 

2 

< Uj^ > , i.e., they both increase at the same rate in the "devel- 
oped" region. 

The one-point velocity-scalar correlations as a function of time 
are shown in Fig. 6-32. As can be seen, there are two non-zero compo- 
nents: 

< Quj^ > < 0U2 > 

___ and -07^— 

This suggests that there are heat/mass fluxes in both the streamwise and 
normal directions. The one-point correlation with the streamwise veloc- 
ity reaches a value of 0.75. The experimental value for this term, 
which was obtained by Tavoularis (1979), is about 0.6. The one-point 
correlation for the normal direction reaches a value of -0.52; the 
corresponding experimental value is -0.4. The difference between our 
computed results and the experimental results for these terms is prob- 
ably due to the Reynolds number difference, as discussed in Section 3A-7 
of this chapter. 

As can be seen from Fig. 6-32, the streamwise velocity and the sca- 
lar quantity are well correlated. This is consistent with the Prandtl 
mixing- length theory, which is based on the assumption that the stream- 
wise velocity component acts as a scalar which is transported verti- 
cally. 

The volume-averaged equations for the scalar intensity and the 
heat/mass flux are given by Eqs. (4-9), (4-11), and (4-12) and are 
discussed in Chapter IV. Figures 6-33 through 6-35 show the time evo- 
lution of each term in these equations. As can be seen, the pressure- 
scalar covariances for both < 0uj > and < 0 U 2 > equations are such 
that they tend to destroy the correlations. This effect is similar to 
that of the analogous terms in the Reynolds shear stress , < U]^U 2 > , 
equation. 
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Figures 6-36 and 6-37 show the 3-D spectra of each term in the 0 
and 0U2 equations. The transfer terms behave in a manner similar to 
that of the transfer term in the equation. 

6.4 Statistics of the Velocity Fields and Models 

As stated above , 24 velocity fields were produced at various Rey- 
nolds and shear numbers. In this section, the effects of these nondl- 
menslonal parameters on various turbulence statistics will be tested and 
their validity will be examined. To do this, we use a power law of the 
nondimensional parameters to fit the data. The accuracy of the fit can 
be judged by r.m.s. relative error obtained from the simulated data and 
the power law. 

A. Reynolds Stress Tensor 

The values of the Reynolds stress tensor components, < >, 

are given in Appendix C. Their behavior is particularly important, 
because they appear in the mean momentum equations and are modeled most 
often. 

The angle of the principal axis of the Reynolds stress tensor rela- 
tive to the computational coordinate system (cf. Eq. (4-17) is given in 
Table 6-5. As can be seen, the principal axis angle varies between 11.5 
and 26 degrees, depending on the shear and Reynolds numbers. The beha- 
vior of this angle as a function of the Reynolds and shear number is 
shown in Figs. 6-38 and 6-39, respectively. As shown in these figures, 
the principal axis angle is Reynolds- and shear-number dependent. It 
decreases with both Re^ and At high Reynolds numbers, however, 

the angle seems to reach an asymptotic value. This value is shear 
number dependent. 

The most common model for < u^Uj > is based on the Boussinesq 
eddy-viscosity concept, Eq. (4-21), i.e., < u^u. > is proportional to 

the mean shear. The problem with this model is that it requires < > 

= < U 2 > = < u^ > for the shear flows. On the other hand, the angle 
of the principal axis of the mean strain, (cf. Eq. (4-19)) is 45 

degrees , while the angle of the principal axis of < > varies 
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between 11.5 and 26 degrees in our flow fields. Thus this model cannot 
represent the Reynolds stress tensor. 

A modified version of this model was suggested by Wilcox and 
Rubesin (1980). Their model has the following form: 


< u.u, 


s - 1 2 ^ 

> = T q <5, 


- S 


2 2 2 

The last term in this model is designed so that < uj » < U 3 > > < up' > 
in a shear flow such as the one we have treated; it makes no 
contribution to < Uj^U 2 > in these flows and this model is identical to 
Eq. (4-21) for this term. For the other components of < u^Uj >, the 
model has the following form: 


. 2 . 1 2 ^ C 

< ^ > = _ q + - 

e 

,2._12_C q^S^ 

< U2 > 2 ^ 2 ' 2 

e 


( 6 - 11 ) 


. 2 . 12 


Wilcox and Rubesin suggested that the coefficient C, for boundary 
layers, is 1.23. From our results, however, the coefficient C is 
found to vary between 0.00022 and 0.0145, depending on the Reynolds 
and shear numbers; the variation of C does not exhibit any clear-cut 
pattern. This suggests that this model is not good for homogeneous 
shear flows. 

A more complicated model for* Reynolds stress tensor was derived by 
Rodi (1976). Rodi's nonlinear algebraic model is based on the assump- 
tion that bj^j is time-independent and is derived for high Reynolds 
number flows. It has the following form: 


< u.u. > 
i J 


T hi 


1 (1-Y) 


- 1 


6 . j ^/e) 


'1 1 + ^ (<^/e - 1) 


( 6 - 12 ) 
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where 


and 


^ - 

ij 


9U. 


au. 

1 


~ ^ u ,U . ^ ^ ~ ^ U .U, ^ -5 — 

i j 3Xj^ j k 3x. 



Cj^ and Y are model constants. Although the values of Y and G]^ 
are constant, this model requires < u| > = < U 3 > in our flow, which 
is not true. 


For the Reynolds shear stress, < Uj^u 2 >, Eq. (4-21) has the fol- 
lowing form: 


3U, 


(6-13) 


where the eddy viscosity, v^, itself needs to be modeled. Models for 
are discussed in Chapter IV. Using the results for each of the 
fourteen simulated flow fields and Eq. (4-23), v,j, = CqL, and Eq. 
(6-13), we calculate the constant C in Eq. (4-23) and correlate it 
with the nondimensional parameters. (Note that the length scale used 
here is the integral length scale, L.) The fitting function used for 
this purpose has the form; 


F = C Re“,^ (6-14) 

where a, 0, and C are determined using least-squares fitting. The 
results are shown in Fig. 6-40. The straight line in this figure rep- 
resents a perfect fit, and the symbols are our simulated results. The 
r.m.s. error shown in the figure indicates the deviation of the data 
from the fit. It is defined as the average of the squares of difference 
between our computed data and the fitting function normalized by the 
average of the data. The solid points are the experimental results of 
Tavoularis (1978). As shown in the figure, v,j,/qL is nearly Reynolds 
number independent (the exponent a in Eq. (6-14) is -0.08) , and is 
approximately proportional to the shear number. 
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(6-15) 



0.265 Re 


-0.08 -1.07 

X 


From this relation, one gets 


< u^U2 > 


2 

q 


0.265 Re 


-0.08-0.07 

X 


(6-16) 


Note that the above relation is good only for > 1. As can be seen, 

O 

< Uj^U 2 >/q is nearly constant in homogeneous shear flows and does not 
depend heavily on Reynolds and shear numbers. 


B. Reynolds Shear Stress Correlation 

< Uj^U2 > 

The values of the Reynolds shear stress correlation, — jp-jp — , 
obtained from the simulated flow fields, are given in Table 6-6. As can 

< Uj^U2 > 

be seen, — , varies from - 0.54 to - 0.74, depending on the 

Reynolds and shear numbers and time. Assuming that it depends only on 

< Uj^U2 > 

Re^ and we use the function (6-14) to correlate p — i — with 

the nondimensional parameters. The results are shown in Fig. 6-41. As 
shown, the correlation decreases slightly with the Reynolds number and 
increases with the shear number. This is expected, as discussed in Sec- 
tion 3A-7 of this chapter. The dependence on Reynolds number is very 
weak. The result is 


“l^2 

u,' u' 

^2 


^ _ -0.043 .0.134 

- 0.515 Re^^ 


(6-17) 


The experimental results of Tavoularis (1979), which obtained at St = 
8 to 11, are also shown in this figure. They lie slightly below our 
computed results. This may be because the Reynolds shear stress corre- 
lation is time-dependent and/or the fitting function is not exact. 

C Reynolds Stress Anisotropy Tensor 

The Reynolds stress anisotropy tensor, j > defined by Eq. 
(4-33). It is a S 3 rmmetric and trace-free tensor, which is zero in iso- 
tropic flows. To determine the effects of Re^^ and ,^s^on the compo- 
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nents of j > again use function (6-14) and least-squares to fit 
the results. The results are shown in Figs. 6-42 through 6-44. As can 
be seen, is almost Independent of Reynolds number. has a 
rather strong effect on bj^j^ and ^ 21 * have much effect 
on bj^ 2 * ^^33 order of magnitude smaller than the other compo- 
nents of b^j and is negligible. It can also be seen from the figures 
that the function does not fit the data well, i.e., there is consider- 
able scatter. It may be because (a) the fitting function is not suit- 
able for j , or (b) b^^ is time-dependent. The latter seems more 
reasonable, since the data figures belonging to one flow fall in a line 
which is not parallel to the fit line. 


D. Pressure-Strain Terms 

As discussed in Chapter IV, the pressure-strain term ‘t’i j > which 
is responsible for transferring energy from one component of the turbu- 
lent kinetic energy to another, can be divided into two distinct parts, 
the fast term and the Rotta term. The magnitudes of each part of the 
pressure-strain tensor are given in Appendix C. Since the fast and 
Rotta terms behave differently, each requires a different model. In 
this section, we test models which were introduced in Chapter IV for 
these terms. 


D-1. The Fast Term 


The part of the pressure strain which is proportional to the mean 
shear is called the fast term (Eq. (4-46)). A model for this term is 
given by Eq. (4-49). This model, which was introduced by Hanjalic and 
Lumley (1974) and W. C. Reynolds (1976), has only one constant Aj^. The 
constant is computed for each component in each flow. Aj^ changes only 
slightly from one flow field to another; the effects of Reynolds and 
shear numbers are very small. However, the constant Aj^ obtained for 
<j>U^ is very different from the constants for the other components of 
. The average value of the constant for each component is shown 
in the first column of Table 6.7. The constant was also fit with the 
function (6.11), and the results are shown in the last three columns of 
this table. The Gibson-Launder model is a simplified version of the 
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above model. It assumes > which is found to be incorrect 
in our flow; therefore this model is poor compared to the more general 
model, Eq. (4-49). 


D-1. The Rotta Term 


The part of the pressure strain which is not explicitly dependent 
on. the mean flow, is the Rotta term (Eq. (4-53)). As shown in 
Table C-5 , the magnitude of 4*^3^ is very small compared to the magni- 
tude of the other components. For this reason, this component is not 
considtjred in this section. 

T»7o models for given in Chapter IV. The Rotta model, 
Eq. (4--55), assumes that the pres sure- strain tensor is proportional to 
the Reynolds stress anisotropy. The constant A 2 for this model is 
calculated using our simulated results and is given in Table 6-8. The 

ij 

function of the form (6.11). The results are given in the table. As 


constant 


A 2 obtained from each component of <j)) 


is fitted by a 


can be seen, 


» varies very much from one flow to another and is not 


well fit as a function of the nondimensional promete rs 
that this model is not likely to be useful. 


This suggests 


Table 6-9 shows the constant A^ and its invariants obtained using 
the modified Rotta model, Eq. (4-56). As shown, A^ has less variation 
and is nearly constant for all of the components. The variation of A^ 
with Rej^ and is obtained using least-square fitting, and the re- 

sults are given in Table 6-9. The constant A^ increases with the 
Reynolds number and decreases with the shear number. If a constant 
value for A^ is to be used, our results suggest A^ = 2.7. 


E« Dissipation Anisotropy Tensor 

The dissipation anisotropy tensor, j , is defined by Eq. (4-57). 
Its magnitude can be calculated from the dissipation tensor, given by 
Table C-3. d^j is a symmetric, trace-free tensor. The dissipation is 
frequently assumed isotropic at high Reynolds numbers, and we expect 
Idj^jl to decrease with Reynolds number. On the other hand, Id^jl is 
zero in isotropic flows when = 0, and one expects it to increase 
with the shear number. We use a least-squares fit to our simulated 
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as a function 


data, using function ( 6 - 11 ) to study the behavior of d^j 
of and Figures 6-45 through 6-47 show the results. Note that 

d 02 very small, so its behavior is not shown here. As can be seen 

from the figures, all of the components decrease with Re^^ and increase 
with This is expected, as discussed above. However, the effect of 

Reynolds number is smaller than expected, suggesting that the dissipa- 
tion may remain anisotropic even at quite high Reynolds numbers. This 
is contrary to the assumption made by many people that the dissipation 
is isotropic at high Reynolds numbers. 

6-5 . Behavior of the Scalar Fields and Their Models 

In this section, the effects of Reynolds, shear, and Prandtl num- 
bers on various turbulence statistics obtained from 56 generated scalar 
fields will be studied. Models introduced in Chapter IV for these terms 
will be tested, and their validity over the range of 
Pr will be examined. The constant(s) for each model will be calcula- 
ted, and their behavior as functions of the non-dimensional parameters 
will be discussed. We use least-squares fits to correlate the results. 

A. The Non-Dimensional Scalar Fluctuations 

The non-dimensional scalar fluctuation, -'*^2 ~ 6 VS 2 L, is a de- 
pendent parameter, as discussed in Section 6-1. Therefore, we suspect 
to be a function of Reynolds, shear, and Prandtl numbers. We use 
the following power-law function to determine its dependence on Re^, 
ydy, and Pr: 

F = C Pr^ (6-18) 

where C, a, 3, and y ate determined by least-squares fitting. The 
results, along with the experimental results, are shown in Fig. 6-48. 

As can be seen, the experimental results behave very similar to the 
simulated results. It can also be seen that ,x ^2 almost independent 
of result is 

^ Pr^-^^^ (6-19) 
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Substituting the definitions of ^<^2 -"^in the above equation, the 

above relation can be approximated as: 


0' q n 0-172 

= __2 pr 

ac/ax^ ^u^/ax^ 


(6-20) 


so the ratio of the r.m.s. scalar fluctuation to the mean scalar gradi- 
ent is proportional to the ratio of, the turbulence velocity to the mean 
velocity gradient. 


B. Turbulence Heat Flux 

The magnitudes of the turbulence heat fluxes , < 0Uj^ > and < 0U2 > , 

obtained from the simulated scalar fields, are given in Appendix C. Note 

that < 0U2 > = 0. As can be seen from Table C-7 , the heat flux ratio, 
< 0Uj^ >/< 0U2 > varies from 1 to 4, depending mainly on the shear num- 
ber. This is expected because, as ,^^<^increases , the flow becomes more 
anisotropic and increases. Therefore, < 0u^ >/< 0U2 > in- 

creases with 


C . One-Point Velocity-Scalar Correlations 

The velocity-scalar correlation is defined as < Ou^^ >/0'u^’ 
(no summation). There are two non-zero components, < 0Uj^ >/0'u^ 
and < 0U2 >/0'u^. They are given in Table 6-10. The magnitudes of 
these terms change only slightly from one flow to another. To correlate 
these terms with the independent parameters of the flow, we use function 
(6-18). The results are shown in Figs. 6-49 and 6-50. They have the 
following form: 


' 6uj^ > 

"FuJ 


= 0.417 Re 


0.021 


,0.138 „ 0.012 
Pr 


( 6 - 21 ) 


and 


< 0Uo > 


n « 0-066 .0.139 ^ -0.061 

- 0.561 Re^ Pr 


( 6 - 22 ) 


Comparing these with Eq. (6-17), we see the one-point velocity-velocity 
and velocity-scalar correlations are very similar. They all increase 
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with the shear number at the same rate and have very weak dependence on 
Reynolds number. 


D. Turbulent Eddy Dlffuslvlty Tensor 
The turbulent eddy dlffuslvlty tensor, 




Is defined by Eq. 


D 


22 


can 


(4-26). It is norr-symmetric. Only two components “^12 and 
be calculated from our 56 simulated scalar fields. The other components 
can be computed when 3C/3x^ and 3C/3x^ are not zero and will be dis- 
Let us consider the behavior of 


cussed later. 


Dj^2 


D 


22 * 


Then we 


study the behavior of the other components using the results obtained 
from two particular runs with non-zero 3C/3x^ and 30/ 3x . Figures 
6-51 and 6-52 show D 22 /qL as power law functions of the 

non-dimensional parameters. The results obtained by using least-squares 
fitting have the following forms; 


*^12 ^ 
qL 


0.259 Re 


-0.014 


-0.72 „ 0.19 
^ Pr 


D 


22 


qL 


„ „ -0.132 -1.1 _ 0.11 

- 0.487 Re, ^ Pr 


(6-23) 

(6-24) 


Using the definition of and and Eq. (6-16), the above rela- 

tions are approximately 

< 6u^ > 

—■ ^ 4— == 0-^78 (6-25) 

< 6u„ > 

— == 0.248 (6-26) 

qo 


These are very similar to the relation for the one-point velocity- 
velocity correlations, Eq. (6-16). This suggests that the eddy viscos- 
ity and dlffuslvlty behave similarly. This is also seen from the exper- 
imental results (Tavoularis (1979)). 


In Chapter IV, we derived a model for the eddy diffusivity tensor. 
The models for Dj ^2 ^22’ have the following forms. 


12 

q^/S 


< > 


Re 


6 P,Y 


(6-27) 
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and 


22 

q^/S 


■ T-ii. i'i 

H'q^ 'V 


+ c. 


Re 


i* . 3' 


Pr 


Y' 


(6-28) 


The magnitudes of the constants in the above relations were ob- 
tained using least-squares of our computed results. The results are 
shown in Figs. 6-53 and 6-54. As shown, the magnitudes of the constants 
are Cj^ = 1.20, C 2 = 0.99, and = 0.569; the constant is 

fixed for both Dj ^2 ^22* constants are nearly Independent of 

Reynolds number and are proportional to Pr^*^^. This model fits the 
data quite well. 

So far, we have discussed only the two components of that 

participate when 30/3x2 is the only non-zero mean scalar gradient. To 
obtain the other components of further cases. Each 

run used the same hydrodynamic field and initial conditions, but one had 
non-zero 3C/3xj^ and the other was for non-zero 30/3x2. Table 6-11 
shows the components of D^j obtained from these simulations. As can 
be seen, has the largest value. This is expected because is 

proportional to Uj^ , which is the largest velocity component. 0^2 and 
D 2 j^ have the smallest values, for similar reasons. Only a few compo- 
nents of are experimentally determined. Table 6-12 shows the 

values of ^ij/^22 obtained from the experimental results of Oorrsin 
(1980) and the present results. The experimental results are obtained 
at Re;i^ = 396 and = 12.76. 

Note that the diffusion equation (2-20) is linear, so the 0 ob- 
tained from different runs can be added or subtracted. Results for a 
case in which the mean scalar gradient is non-zero in all three direc- 
tions can be obtained from three runs , each of which has a mean gradient 
in one direction. However, the 8 fields are not uncorrelated, because 
the same velocity field is the same. Then, < ^ 1^2 ^ 
necessarily zero when and 62 are results from two different runs. 


E. Turbulent Prandtl Number 

The turbulent Prandtl number, 
to eddy diffuslvity (Eq. (4-36)). 


Pr^, is the ratio of eddy viscosity 
The magnitude of turbulent Prandtl 
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number for each flow field is given in Table 6-13. Its behavior as a 

function of Reynolds, shear, and Prandtl numbers is shown in Fig. 6-55. 

As can be seen, Pr^ is almost Independent of Re^ and It de- 

creases when molecular Prandtl number Increases, and vice versa. The 
result of least-squares fitting of the data is: 

Pr^ == 0.528 pj.-0.105 (6-29) 

TKe^ v'aIol^ of 

Pr^ obtained from 56 simulated flow fields were used to test models 
described in Chapter IV. The results are shown in Figs. 6-56 through 
6-59. 

The first three of these figures test Reynolds', Wassel's, and 

Crawford's models, given by Eqs. (4-37) through (4-39), respectively. 
All of the models behave similarly and fit the simulated results rather 
poorly. The suggested values for the model constants suggested in the 
literature are different from the constants obtained from our simulated 
results. 

Figure 6-59 tests the model derived in Chapter IV, Eq. (4-40). 
This model fits the data much better than the models mentioned above. 
The model can be approximately simplified to: 

Pr^ = 2 

^0.2 pj.0.12 ^j^^20 ^ — + 0.169^ 

F. Pressure-Scalar Covariance 

Like the pressure-strain, the pressure-scalar covariance contains 
two parts. The magnitudes of each part of the pressure-scalar covari- 
ance computed from our simulated flow fields are given in Appendix C. 
In this section we examine the pressure-scalar covariance models intro- 
duced in Chapter IV. 
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F-1. The Fast Term 


The part of the pressure-scalar covariance which is proportional to 
the mean velocity gradient called the fast or rapid term 

(Eq. (4-60)). 

Let us consider the model given by Eqs. (4-69) and (4-70). The 
constant in this model is 2/5, according to the assumptions made in 

Chapter IV. In this section we use our simulated results and calculate 
A^. The behavior of A^ as a function of the nondimensional parameters 
is obtained using the least-squares fitting and function (6-21). The 

results are shown in Jigo t — 6- ' 70 to — 6 " 71-“ a nd — a - l - se * in Table 6-14. As 

shown, the value of A^ changes very slightly throughout the 56 flow 
fields (i.e., the variance of A^ is small) and has small dependency 
on Re;)^, and Pr. This suggests that much of the important physics 

is captured by the model. On the other hand, the computed value of A^ 
is less than one-half of the expected value. This suggests that the as- 
sumptions (mentioned in Chapter IV) made in deriving the model are not 
sufficient, and a more complicated model may be needed to get better 

results. 

Let us now consider Luinley's model, Eqs. (4-72) and (4-73). This 

model satisfies all of the conditions of the previous one plus a new 

condition mentioned in Chapter IV. The constant A^ in this model is 
predicted to be 1/2. The values of A^ and its behavior as a function 
of Re^, and Pr for the two components of 4*^^^ * given in 

Table 6-15. As shown, the variation of A over 56 flow fields is 
small, and therefore the model captures much of the important physics. 
However, the magnitude of the constant is different from 1/2 and varies 
from one component to another. 

F-2 . The Rotta Term 

The second part of the pressure-scalar covariance, is 

defined by Eq. (4-74). This part does not explicitly depend on the mean 
flow. Equations (4-79) and (4-80) represent a model for this term. The 
model is analogous to the modified Rotta model used for the pressure- 
strain covariance. The constant Aj^ obtained by using our simulated 

results is given in Table 6-16. The behavior of Aj^ as a function of 
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Re^ and are also given in the table. As we see, the constants for 
both components of are nearly the same and have very small vari- 
ation due to the and Pr. This is similar behavior to that 
of the modified Rotta model for the pressure-strain terms. 

6.5 Comparison with the Other Flows 

Nearly all turbulent flows in nature are inhomogeneous, i.e., the 
shear varies through the flow. The homogeneous flows may be structur- 
ally different from their inhomogeneous counterparts; the two kinds of 
flows might require different types of models , and thus our results and 
discussion for mixing of a passive scalar in homogeneous flows might not 
be valid for inhomogeneous flows. We shall briefly compare some of our 
results with those for inhomogeneous flows. 

Table 6-17 compares non-dimensional quantities obtained in the log- 
arithmic region of a boundary layer and in one of our runs. These 
values are nearly constant in the logarithmic region but differ in the 
outer region. 

The shear number was found to play an important role in homogeneous 
flows. Its value as a function of y^ in the boundary layer is given 
in Table 6-18. It is nearly constant in the logarithmic region but 
changes rapidly in the other regions. In this region, the agreement 
with our results is quite good. In the inner region, the shear number 
changes rapidly and the wall has a very strong effect on the behavior of 
the turbulence. Both of these effects make the flow very different from 
a homogeneous flow, and our results cannot be applied. Clearly, spe- 
cific models which account for these effects are needed and are, of 
course, in common use. In the outer (wake) region of the boundary 
layer, the shear number rapidly becomes small and there are again 
inhomogeneity effects which make comparison with our results very 
difficult. The interaction with the irrotational outer flow also has a 
strong effect. 

Table 6-19 shows the variation of some of the turbulence statistics 
through a turbulent free shear layer and compares them with our homoge- 
neous shear flow results. The turbulence statistics are nearly constant 
for -0.5 < n < 0.5. In this region the shear rate, 8Uj^/3x2, is 
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nearly constant and the results, as expected, are comparable to our 
homogeneous shear flow results. The shear number is also nearly con- 
stant in this region and is of the same order of magnitude as in bound- 
ary layers. In the region where ln| > 0.5, the interaction with the 
irrotattonal flow causes strong inhomogeneity effects and the comparison 
with our flow is difficult. 
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Table 6-1 


VELOCITY FIELD 


Parameter 

Notation 

Dimension 

Turbulence Velocity 

q = < u^Uj: 

L T"^ 

Kinematic Viscosity 

V 

l2t-1 

Integral Length Scale 

L 

L 

Taylor Microscale 

X 

L 

Time 

t 

T 


Table 6-2 
SCALAR FIELD 


Parameter 

Notation 

Dimension 

Turbulent Scalar Intensity 

6' = < 

0 

Diffusion Coefficient 

D 

L^T“^ 

Mean Scalar Gradient 

3C’ 

^ -1 

^2 ■ 



0 L ^ 

Integral Length Scale 

Le 

L 

Taylor Microscale 


L 
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TABLE 6-3 

Ransie of Non-Dimensional Parameters 


[parameters 

Present Data 

Experimental Data 


10 - 120 

130 - 400 


5-70 

3 - 14.8 

Pr 

0.2 - 5 

0.7 

St 

0-7 

0-12 


TABLE 6-4 

Non-Dimensional Parameters Por Each Flow Field 


Simulations 

Parameters 

St=4 


St = 6 

St = 7 


Re 

44.560 

52.057 

59.215 

66.762 

HS64A 


18.549 

18.983 

19.485 

19.993 


Pr 

1.000 

1.000 

1.000 

1.000 


Re 

18.996 

21.840 

24.545 

27.175 

HS64B 


23.727 

25,178 

26.790 

29.397 


Pr 

1.000 

1.000 

1.000 

1.000 


Re 

60.348 

75.711 

92.476 

10.080 

HS64G 


66.879 

69.205 

68.924 

68.878 


Pr 

1.000 

1.000 

1.000 

1.000 


Re 

24.202 

27.626 

31.233 

35.222 

HS64J 


11.228 

12.617 

14.493 

16.328 


Pr 

1.000 

1,000 

1.000 

1.000 


Re 

29.156 

34.132 

39.054 

44.550 

HS64M 


8,936 

10.333 

11.987 

13.853 


Pr 

1.000 

1.000 

1.000 

1.000 


Re 

64.310 

71.840 

78.990 

86.910 

HS640 


6.376 

7.120 

7.738 

8.236 


Pr 

1.000 

1.000 

1.000 

1.000 


Re 

18.996 

21.840 

24.545 

27.175 

HS64C 


23.727 

25.178 

26.790 

29.397 


Pr 

5.000 

5.000 

5.000 

5.000 
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TABLE 6-4 


Continued 


Simulations 

PArameters 

St=4 

St=5 

St=6 

St=7 


Re 

44,560 

52.057 

59,215 

66,762 

HS64D 


18,549 

18.983 

19,485 

19.993 


Pr 

5,000 

5,000 

5.000 

5.000 


Re 

44,560 

52,057 

59.215 

66,762 

HS64E 


18,549 

18.983 

19.485 

19,993 


Pr 

0,200 

0,200 

0,200 

0.200 

HS64F 

Re 

18,996 

21.840 

24.545 

27,175 


23,727 

25,178 

26.790 

29.397 


Pr 

0,200 

0.200 

0.200 

0,200 

HS64H 

Re 

60,348 

75.711 

92,476 

10,080 


66.879 

69,205 

68.924 

68.878 


Pr 

5.000 

5.000 

5.000 

5,000 


Re 

60,348 

75,711 

92,476 

10,080 

HS64I 


66,879 

69,205 

68,924 

68,878 


Pr 

0,200 

0,200 

0,200 

0,200 


Re 

29,156 

34.132 

39,054 

44,550 

HS64N 


8.936 

10.333 

11.987 

13,853 


Pr 

0.200 

0,200 

0,200 

0.200 


Re 

64,310 

71,840 

78.990 

86,910 

HS64P 


6,376 

7.120 

7,738 

8,236 


Pr 

0,200 

0,200 

0.200 

0.200 
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TABLE 6-5 


Principal Axis An^le ( Radians) 
Eauation 4-17 


Siniiul atlon 

ST = 4 

ST=5 

.SI.=.i^ 

ST=7 _ 

HS64A 

0.31886 

0.27182 

0.23674 

0.22285 

HS64B 

0.32776 

0.27180 

0.22677 

0.20689 

HS646 

0.26261 

0.21147 

0.17420 

0.15333 

HS64J 

0.40883 

0.34481 

0.29152 

0.26612 

HS64M 

0.43080 

0.37272 

0.32677 

0.30371 

HS640 

0.45910 

0.42296 

0.40323 

0.39940 


TABLE 6-6 


One-Point Correlations of the Velocity Field 




UJ 


u 


2 


Simulation 

St = 4 

St=5 

St = 6 

St=7 

HS64A 

0.68350 

0.67890 

0.68780 

0.67250 

HS64B 

0.67680 

0.67560 

0.70390 

0.69340 

HS64G 

0.72830 

0.73470 

0.76240 

0.74390 

HS64J 

0.62960 

0.63344 

0.63720 

0.61760 

HS64M 

0.61800 

0.61 666 

0.60490 

0.59520 

HS640 

0.54260 

0.54720 

0.54800 

0.54280 
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TABLE 6-7 

Coristant For the Pressure-Strsin Models Eq. 4-49 

The Fast Term 


GL 6 

Fitting Function = C Re 


E Q u a t i o n 

Ai 

C 

a 

e 

r.ni.s. error 

MH 

W 

-6.94 ±0.519 

9.75 

-0.028 

-0.082 

0.038 

Mj 

J = 2 

-1.83 ±0.400 

1.55 

0.038 

0.008 

0.060 

i = 3» 

J = 3 

-1.57 ±0.380 

6.82 

-0.185 

-0.282 

0.093 

i = l > 

J = 2 

-1.52 ±0.077 

2.01 

-0.082 

0.010 

0.032 


TABLE 6-8 

Constant For the Pressure-strain model» Eg. (4-53) 

Rotta Model 

FittinS Function Re*^ 


Eauation 

A 

C 

a 

fi 

r.iti.s. error 


0.673 ±0.599 

1 .143 

0.678 

-1.195 

0.15 


0.673 ±0.570 

0.734 

0.715 

1)0 

c- 

0 

1 

0.13 

tlBM 

0.679 ±0.452 

0.238 

0.813 

-0.776 

0.15 


TABLE 6-9 

Constant For the Pressure-Strain Models Eq. (4-56) 

Modified Rotta Model 


FittinS Function A =C Re^ ^ ^ 

3 


E Q u a t i o n 

A j 

C 

a 

6 

r.m.s. error 

i=l» J=1 

2.862 ±0.349 

3.342 

0.144 

-0.225 

0.045 

i=2» J=2 

3.080 ±0.428 

4.200 

0.119 

-0.245 

0.048 

i=^l» J = 2 

2.472 ± 0.426 

1 . 633 

0.330 

-0.266 

0.085 
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TABLE 6-10 

One-Point Veloci ta-Scalar Correlation 


<u^ 0 >/ 6’ 


SJjnulation 

I 

St = 4 

St = 5 

Bt=6 

St=7 

HS64A 

1 

0,667 A 

0.6513 

0.6351 

0.5822 

2 

0.7064 

0.6742 

0.6370 

0.5891 

HS64B 

1 

0.6504 

0.6130 

0,5726 

0.5197 


2 

0.7476 

0.7302 

0.7072 

0.6410 

HS64G 

1 

0.7949 

0.7541 

0.7041 

0.6583 


2 

0.7803 

0.7746 

0.7765 

0.7320 

HS64J 

1 

0.8552 

0.8236 

0.7709 

0.7102 


2 

0.6890 

0.6542 

0.5878 

0.5197 

HS64M 

1 

0.7721 

0,7662 

0.7687 

0.7244 

2 

0.6570 

0.6134 

0.5570 

0.5222 

HS640 

1 

0.8131 

0.8058 

0.8031 

0.7564 

2 

0.5682 

0.5506 

0.5195 

0.4999 

HS64C 

1 

0.7877 

0.7464 

0.6871 

0.6524 


2 

0.6674 

0.6513 

0.6351 

0.5822 

HS64D 

1 

0,7090 

0.6724 

0.6291 

0.6078 


2 

0.6504 

0.6130 

0.5726 

0.5197 

HS64E 

1 

0.6561 

0.6806 

0.7051 

0,7414 


2 

0.7949 

0.7541 

0.7041 

0.6583 

HS64F 

1 

0.6054 

0.6129 

0.6282 

0.6706 

2 

0.8552 

0.8236 

0.7709 

0.7102 

HS64H 

1 

0.7315 

0.7887 

0.8301 

0.8614 

0 

0.7721 

0.7662 

0.7687 

0.7244 

HS64I 

1 

0,7306 

0,7770 

0.8099 

0.8366 

2 

0.8131 

0.8058 

0.8031 

0.7564 

HS64N 

1 

0.5949 

0.6310 

0.6506 

0.6805 

2 

0.7877 

0.7464 

0.6871 

0.6524 

HS64P 

1 

0.5495 

0.5793 

0.5991 

0.6198 

2 

0.7090 

0.6724 

0.6291 

0.6078 
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TABLE 6-11 

Components of Nomelized Eddu Diffusivity Tensor 


Eciuat ion 

St=3 

St=4 

St«5 

St=6 

i = l r 

j=l 

-0.08770 

-0.09750 

-0.10690 

-0.09810 

i=2r 

J=2 

-0.03200 

-0.02530 

-0.01792 

-0.01200 

i=3» 

J=3 

-0.03430 

-0.03325 

-0.03335 

-0.02590 

i = l » 

J=2 

0.04645 

0.04786 

0.04696 

0.03910 

i=2» 

j=l 

0.03253 

0.03026 

0.02462 

0.01910 


TABLE 6-12 

Experimental and Computed Malues of Normalized 
Turbulent Eddy Diffusivity 


Eouations 

Computed Oalues 
at St=3 

Experimental Results 
(Corrsin (1980)) 

i=lf j=l 

-6.00 

* 

1=2» j = l 

1.40 

* 

i=l» j=2 

2.30 

2.20 

i=2» j=2 

1.00 

1.00 

i=3» j=3 

-1.60 

-1.60 


t The experimental values are not available* 
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TABLE 6-13 


Turbulent Prandtl Number 
Eeuation <4-36> 


Simulations 

St=4 

St=5 

St=6 

St=7 

HS64A 

0.6617 

0.7089 

0.7850 

0.8585 

HS64B 

0.5984 

0.6173 

0.6766 

0.7550 

HS64G 

0.6589 

0.6950 

0.7445 

0.7945 

HS64J 

0.5800 

0.6125 

0.7006 

0.8010 

HS64M 

0.6093 

0.6602 

0.7321 

0.7950 

HS640 

0.6143 

0.6386 

0.6917 

0.7269 

HS64C 

0.5302 

0.5452 

0.5900 

0.6522 

HS64D 

0.6129 

0.6538 

0.7192 

0.7870 

HS64E 

0.8609 

0.9190 

1.0280 

1.1200 

HS64F 

0.8769 

0.8910 

0.9996 

1.1310 

HS64H 

0.6461 

0.6790 

0.7241 

0.7695 

HS64I 

0.7217 

0.7716 

0.8417 

0.9125 

HS64N 

0.8584 

0.8722 

0.9208 

0.9710 

HS64P 

0.7482 

0.7712 

0.8166 

0.8335 
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TBALE 6-14 


Constant for Pressure-Scalar Covariance Models Eqs» (4-69) and (4-70) 

The Fast Term 


Fittind Function = c Re^ Pr 


E « u a t i o n 

^0 

C 

a 

3 

Y 

r.m.s. error 

(4-69) 

0.229 +0.079 

0.543 

0.025 

-0.341 

-0.096 

0.17 

(4-70) 

0.144 +0.066 

0.312 

-0.018 

-0.287 

-0.148 

0.11 


TABLE 6-15 

Constant for Pressure-Scalar Covariance Model? Eq«( 4-71) 
Lumley's Model for the Fast Term 


Fittind function Aq = C Re“ ^ Pr ^ 


i 

^ 

C 

a 

3 

Y 

r.m.s. error 

1 

0.676 +0.081 

0.754 

0.014 

-0.059 

-0.043 

0.17 

2 

0.276 ±0.103 

0.085 

-0.035 

o 

o 

-0.018 

0.18 


TABLE 6-16 

Coristatnt for the Pressure-Scalar Model 
The Rotta-Like Term 


Fittind Function Aj^ = C Re°^ Pr ^ 


E a u a t i o n 



C 

a 

8 

Y 

r.m.s. error 

(4-79) 

-1 . 580 + 0 . 677 

1.031 

0.170 

-0 . 128 

-0.307 

0.16 

(4-80) 

-1 .789 ±0.730 

0.744 

0.461 

-0.333 

-0.157 

0.14 
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Table 6-17 


Comparison of Turbulence Statistics In Homogeneous Shear Flows 

and Boundary Layers 


Quantity 

Homogeneous Shear Flow 
Run HS64J 
at St = 7 

4? 

Boundary Layer 

< > 

2 

q 

-0.147 

-0.130 

< Uj^e > 
q0' 

0.547 

0.50 

< u^e > 

q0' 

-0.160 

-0.20 

Pr^ 

0.80 

1.05 

^12 

“22 

-2.30 

-2.40 


*Johnson (1959) 
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Table 6-18 


Variation of Non-Dimensional Quantities in Channel Flow^ 

«** 
u 6 

Re. = -£ = 13800 

0 V 


Sublayer 


>Log-region 
Centerline 

%oin (1981). 

Rej^ is Reynolds number based on the 
Integral length scale. 

6 is channel half-width. 

u^ is mean velocity at the centerline. 
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Table 6-19 


Comparison of Turbulence Statistics in Homogeneous Shear Flow 

and Free Shear Layers 


Quantity 

4r 

Free Shear Layer 

Homogeneous Shear 
Flow , Run HS64 J 
at St = 7 

n = -1.0 

n = -0.75 

n = 0.5 

n = 0.8 

n = 0.5 

n = 0.75 

n = 1.0 

n 

0.312 

0.38 

3.45 

0.50 

r\ /. n 
u •♦tu 

0.28 

n 1 

0.62 

< Qu^ > 
0'u^ 

0.55 

0.59 

0.60 

0.65 

0.42 

0.332 

0.23 

0.520 

Pr^ 

0.60 

0.62 

0.65 

0.82 

0.79 

0.70 

0.56 

0.80 


0.987 

1.16 

1.25 

1.58 

2.00 

1.82 

1.82 

16.00 


6231 

8107 

9340 

1 

9262 

j 

8936 

7678 

- - - - - - 1 

6234 

140 


*Batt (1977). 

The values obtained at — Xj^q = 18.5" for = 23 ft/ sec and U 2 = 0.0. 


^ 12(x2 - 


^1 " *10 


, Re, 


UiL 
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where L is integral length scale 
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Fig. 6-1. Nondlmensional scalar quantity, ^g2> ® function of tine for 

various runs with different Initial values. 



Fig. 6-<:. Integral length scale ratio of the velocity and scalar field as a 
function of time for various runs with different Initial values. 
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Fig. 6-4. Three-dimensional spectra of the turbulent kinetic energy compo- 
nents. 
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Fly. 6 - 5 . Turbulent kinetic energy and Its components as a function of title. 



Fly. 6-6. Integral length scales in the »^-direction as a function of time. 
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F(g. 6-7. Taylor microscales in the Xj -direction as a function of time. 
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‘icf id 

WAVE NO. 

Fig. 6-li;. One-dimensional spectra of velocity components and scalar in the 
xpdirection. 
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Fi^* 6-13. One-dimensional spectra of velocity components and scalar in the 
x^-direction. 



Fiy. 6-14. One-dimensional spectra of velocity components and scalar in the 
Xj-direction. 
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FiQ* Three-dimensional spectra of two conponents of pressure and total 

pressure at St » O.U, 



ICf id 


WAVE NO. 

Fig. ti-16. Three-dimensional spectra of the components of pressure and total 
pressure at St = 4.0. 
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Fly. 6-17. One*po1nt velocity correlations as a function of time. 
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Fiy. b-18. Components of Reynolds stress anisotropy tensor as a function of 
time. 
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Fly. b-19. Various terms in the < u( > equation as a function of time. 
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St 

Fig. 6~Z0. Various terns in the < > equation as a function of time. 
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Fia. Various terms in the < uj > equation as a function of time. 
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Fig. 6-25. Three-dimensional spectra of various terras In the equation as 
a function of time. 



Fig. 6-26. Three-dimensional spectra of various terms In the u| equation as 
d function of time* 
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Fla. S-27. Three-dimensional spectra of various tenns In the UiUj, equation 
as a function of time. 



Fly. 6-28. Three-dimensional spectra of various tenns In the u^u^/2 
equation as a function of time. 
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Fig. 6-30. Three-dimensional spectra of the scalar fluctuations at different 
times. 
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Fig* 6-33. Various terms in the < 0^ > equation as a function of time. 
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Fig, 6-34* Various terms in the < 6u;^ > equation as a function of time* 
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Fiy. 6-37. Three-dimensional spectra of various terms in the 6uj; equation. 
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F1g. 6-38. The angle of the principal axis of the Reynolds stress 
tensor as a function of Reynolds number ( Re^) 
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Fig. 6-39. The angle of the principal axis of the Reynolds 
stress tensor as a function of shear number 
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Fiy. 6-4U. The constant for the eddy viscosity model (C vj/qL) as a 
function of the fittiny function. 



Fiy. 6-41. Reynolds shear stress correlation as a function of the fittiny 
function. 


122 






123 







124 





020 0.21 024 OK 028 0 20 022 024 0.28 028 0.40 0.42 044 0.48 


Fij. 6-46. Dissipation anisotropy, d^2« * function of the fittiny 

function. 
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Fiy. b-bZ. Normalized turbulent eddy diffusivity, D^^/qL, as a function of 
the fitting function. 
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Fig. 6-53. Turbulent eddy dlffuslvity, *5 « function of its model. 
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Fig. 6-b9. Turbulent Prandtl number as a function of the new model. 
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Chapter VII 

CONCLUSIONS AND RECOMMENDATIONS 


7 .1 Conclusions 


In this work we have numerically solved the three-dimensional , 
time-dependent equations of motion and the diffusion equation for two 
cases , homogeneous isotropic turbulence and homogeneous sheared turbu- 
lence . In order to do this , we were required to work at rather small 
Reynolds numbers. The results obtained have been compared with avail- 
able experimental results, and good agreement was obtained. 

For the homogeneous isotropic flows , it has been found that the 
Reynolds and Prandtl numbers (within the range considered) have strong 
effects on the decay rates and the length scale ratios. The decay expo- 
nent of the turbulent kinetic energy, n, decreases with Increasing 
Reynolds number (from 2.5 at Re^ = 0 to 1.2 at Re^ ~ 25). The decay 
exponent of the scalar fluctuation intensity, m, decreases with both 
Reynolds and Prandtl numbers. Both the decay exponent ratio, R, and 
the length scale ratio, Increase with the Reynolds number when 
Pr < 1 and decrease with the Reynolds number when Pr > 1 . The abso- 
lute value of the velocity derivative skewness, Sk, has a maximum 
value of 0.5 at Re^ “ 20. Sk is zero at Re;^ = 0 and reaches an 
asymptotic value of approximately 0 .4 as Rej^ approaches 100 . 

In homogeneous shear flows , unlike homogeneous Isotropic flows , 
the initial conditions, such as the initial length scale ratio, ^u/^ 0 > 
and the initial non-dimensional scalar fluctuation do not affect the 
scalar field. Many of the turbulence statistics are nearly independent 
of the Reynolds number. The shear number has the most important influ- 
ence in this flow. In homogeneous shear flows, the length scales grow 
faster than in the Isotropic flows. The behavior of the dissipation (or 
transfer) length scale and integral length scales are similar, but the 
mixing length seems to behave differently. The turbulent kinetic energy 
first decreases and then increases. Due to the shear, the flow is ani- 
sotropic; the streamwise velocity fluctuations are the largest, while 
the vertical velocity fluctuations are the smallest. The dissipation 
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tensor is very anisotropic and does not become isotropic even at rather 
high Reynolds numbers. This is contrary to the assumptions made by many 
people. The behavior of the scalar fluctuation intensity is very simi- 
lar to that of the streamwise turbulent kinetic energy components; they 
are well correlated, and their spectra are similar. It is unclear 
whether this flow reaches statistical equilibrium. Our results show 
that < Uj^U2 >/q , < 0uj^ >/q0'» and < 0U2 >/q0' are nearly constant 
in this flow. However, the other components of bj^j slightly change 
with time. 

From the three-dimensional spectra of the production, dissipation, 
pressure-strain, and the transfer terms, it is found that influence of 
the transfer term is very small compared to that of the other terms . In 
sheared turbulence, this term, contrary to a common assumption, trans- 
fers energy from both large and small eddies to the moderate-sized 
eddies . 

The three-dimensional spectra of the components of pressure field 
(p^^^ and p^^^) have been studied. It has been found that the part of 
pressure which is explicitly associated with the mean flow dominates at 
the large scales and has a k~^ shape at moderate wave numbers . The 
second part of the pressure is important at the small scales . The 
relative magnitude of the two parts of the pressure depends on the shear 
rate . 

The turbulent Prandtl number is almost independent of the Reynolds 
number, and decreases as the molecular Prandtl number increases. A new 
model for the turbulent Prandtl number seems to correlate the data bet- 
ter than a number of other models . 

Several models for various turbulence statistics were tested, and 
the following conclusions made : 

The eddy viscosity model fits the turbulent shear stress quite 
well, but is incapable of predicting the diagonal components of the 
Reynolds stress accurately. The parameter in the model is proportional 
to the inverse of the shear number. Rodi's algebraic model for the 
Reynolds stress is also incapable of predicting the diagonal components 
of the Reynolds stress in this flow, but its coefficient is relatively 
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constant for the turbulent shear stress. Eddy diffusivity models behave 
similarly to the eddy viscosity model. Their parameters are nearly pro- 
portional to the Inverse of shear number and depend slightly on Prandtl 
number. 

The Rotta model for the pressure-strain term correlates the data 
rather well. The modified Rotta model correlates the data even better. 
The constant for this model is about 2.7. Similarly, the model similar 
to the modified Rotta model for the pressure-scalar covariance fits the 
data very well. The models for the fast term are not good and need to 
be improved . 

7 .2 Recommendations 

In this work, we have assumed the scalar is passive; therefore, we 
have neglected the effects of buoyancy. It would be interesting to 
study the effects of buoyancy in these flows . Our simulated flow fields 
have provided "exact" results for the flow fields. We have used these 
results to check some of the models for various volume- averaged turbu- 
lence statistics. We suggest using the "exact" results to obtain some 
models for the subgrid scales of the scalar fields . 

Finally, it should be possible to extend this work to simulate the 
reacting flows. , 



Appendix A 


LINEAR STABILITY ANALYSIS 


In this section we study the stability and accuracy of numerical 
methods of solving the equations of motion (2-18) and the diffusion 
equation (2-20) . 

For the sake of simplicity, let us assume that the problem is two- 
dimensional. In two dimensions, the Fourier transform of the linearized 
equation (2-18) has the form: 


+ v(-k.k. + 2St k,k„ - S^t^k?) u- 
j j 12 r Z 


(A-1) 


where u ^ is assxmied to be constant. In component form, this becomes 
3ui " r 

^ = u^ |^(- u^jk^ . “o2‘'l^ ‘ ‘'1"2 


3t 


+ St u^^k, ) i + v(-k^k^ + 2 St k, k, - S^t^kJ) 


- ikj^p - S u^ 


(A-2) 


au^ 

9t 


= u„ (- u .k, 

2 L J 


+ St u „k, ) 1 + v(-k.k. + 2 St k-k- - S^t^k 
o2 1'' J J 12 


J>] 


- ik2p + St ikj^p 


(A-3) 


Let us define 


a = (- ^ ^o2*^l^ ^ ^ ~ ^2 ^ ^ *^1^2 ~ 

(A-4) 

which is simply a known function of the wavenimiber. Then the equations 
(A-2) and (A-3) become 
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9 A A A 

IT “ 

A 

^ 

= u^a - Ik^p + St k^p 


(A-5) 


Following the standard linear stability analysis for Navier-Stokes 
equations , since the dynamics of the pressure is totally dependent on 
the dynamics of the velocity field, we consider the following equations. 



auj^ + Su2 


(A-6) 



(A-7) 


The solution of the above equations using the Adams-Bashforth method 
(2-11) is 

“'n+l ''n , 1 ../Q ^'n "‘n-1. , 1 ... /a o\ 

U2 ~ “2 2* ~ ®“2 ^ 2" ~ ^^^2 ^ (A~8) 

Let us assume 


= Ru" (A-9) 

First let us consider the stability of the U 2 equation. From (A-9) 
and (A-8), for 1=2, one can get: 

R = l+|-At3a-y (A-10) 


or 

- (1 + j Ata) R + a = 0 (A-11) 

or 

„ 1 3At . .n 9 .X 2 ~ 

where a is defined in (A-4). 
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We can write the maximum value of a as 


a = (iN(, + Nj))/At 

where is the Courant number and is defined as 


= (lu.l k. + St|u„| k, ) At 
C j max j max 2 max 1 max 


and Ny is the viscous stability parameter 


(A-13) 


(A-14) 


N 


^ + 2 St k, k„ - S^t^k^ 

j max I max 2 max 


= f-k. k 

V L J 

For our problem the maximum wavenumber is 


2 

1 max_ 


vAt (A-15) 


max ^2 max ~ % max ~ ^max 


where A is mesh size. Note that the computational box length is 2ir , 
so A = 2ir/(N-l), where N is the number of mesh points in each direc- 


tion. Therefore, k „ becomes: 
’ max 


max 


N-1 


(A-16) 


From (A-14), (A-15), and (A-16), the values of N and Ny are 


N, 


V" + St) lu, 1 1 At (A-17) 


and 


■ (^) [■ 


3 + 2 St - S^t^ I vAt 


2 


(A- 18) 


and are time dependent. (Note that St varies between - V 2 and V 2 ) • 

Following the same procedure , one can set the following N^. and 
Ny for the Uj^ equation: 

(3 - 2 St - + y S ^ 


vAt 


max 


(A-19) 
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and is the same as before, Eq. (A-14). The magnitudes of and 

R 2 of Eq. (A-12) depend on the magnitudes of N^, and N^. Rj^ is 
the physical root and R 2 is the spurious root. It should be noted that 
the method is stable even if |r| > 1, because in shear flows the vel- 
ocities increase with time. The limit on R is obtained by comparing 
Eq. (A-9) and the analytical solution of Eqs. (A-6) and (A-7) . 

The magnitudes of Rj^ and R 2 for various and (Eqs . 

(A-12) and (A-13)) are given in Table 2-1. 

For the diffusion equation (2-20), we follow a similar procedure, 
and the Important parameter governing stability is : 

Njj = D (- 3 + 2 St - S^t^) At (A-20) 

For the three-dimensional case, the Courant number is defined as 

Nr = +(1+St)|u, I + Iu_| I At (A-21) 

C 2 |_ 1 max 2 max 3 maxj 

Stability is affected by St and S|u„/u, 1 , but as we keep 

2 1 max 

I St I < 0.5, the effect is small. Also, as we increase S, lu 2 /uj^| 

becomes small and therefore S|u„/u, | is always small compared to 

2 1 max 

the other term in Eq . (A-19) . 

The time step at which the results are reasonably accurate and sta- 
ble is chosen by using Table 2-1 . We chose N^. = 0 .1 to 0 .3 and Ny = 
0.1. Note that N^ = N^/Pr. 


139 



Appendix B 

COORDINATE TRANSFORMATION 


In this section, we consider the coordinate transformation which 
was used to eliminate the non-periodic terms in the equations of motion 
and the diffusion equation. 

The transformed coordinates are 


xj^ = - SX2T 

X£ = X2 

X3 = X3 

t = T 


(B-1) 


where X^ and T are the physical (laboratory) variables. As shown 
below, the tranformed coordinate moves with time. 



Using the chain rule, the derivative of any function with respect 
to X^ and T can be obtained as follows: 



(B-2) 

(B-3) 

(B-4) 
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d 

3T 


at 


- Sx 


2 3x, 


3x5 


3 x„ 


- 2 St 


ax. 


,2 

3 X, 


3x, 3x, 

3 


^‘--2 
2 


^^2 2 
+ S t 


,2 

3 X, 


3x, 


(B-5) 

(B-6) 

(B-7) 

(B-8) 


The above relations were used to transform the non-periodic 
equations (Section 2.5b) into Eqs. (2-17) through (2-20). 
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Appendix C 

TABUIATED DATA FOR THE SIMULATED FLOW FIELDS 


In this appendix, we present a set of tables containing various 
information, raw data, and turbulence statistics obtained from 24 
simulated homogeneous shear flow fields and 56 scalar fields. 

The Initial values used In all of the simulations are: 

Turbulence velocity, q = 0.600. 

Integral length scale, L = 0.71. 

Taylor micro-scale, X = 0.22. 

r.m.s. scalar fluctuations, 0’ == 0.600. 

Integral length scale of the scalar field, Lq = 0.51. 

Taylor microscale of the scalar field, Xq = 0.17. 

The above parameters are nondimenslonallzed , as stated In Section 

2.3. 

The following table shows other parameters used for each simula- 
tion. 
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TABLE C-1 

Dimensional Parameters For Each Simulation 


Simulation 

Shear Rate 
S 

Viscosits 

V 

Dif f usi vi ts 
D 

Mean Scalar 
Gradient 

HS64A 

8.851 

0.00500 

0.00500 

8.500 

HS64B 

8.851 

0.01000 

0.01000 

8.500 

HS646 

42.890 

0.00500 

0.00500 

42.500 

HS64J 

2.860 

0.00525 

0.00525 

2.860 

HS64M 

1.716 

0.00350 

0.00350 

1.716 

HS640 

2.230 

0.00175 

0.00175 

2.230 

HS64C 

8.851 

0.01000 

0.00200 

8.500 

HS64D 

8.851 

0.00500 

0.00100 

8.500 

HS64E 

8.851 

0.00500 

0.02500 

8.500 

HS64F 

8.851 

0.01000 

0.05000 

8.500 

HS64H 

42.890 

0.00500 

0.00100 

42.500 

HS64I 

42.890 

0.00500 

0.02500 

42.500 

HS64N 

1.716 

0.00350 

0.01750 

1.716 

HS64P 

2.230 

0.00175 

0.00875 

2.230 
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TABLE C-2 


Components of the Reynolds Stress Tensor 


Simulation 

Eauat ion 

St = 4 

St=5 

St=6 

St=7 


i=l. 

j =1 

0.27469 

0.32620 

0.37760 

0.42920 

HS64A 

i=2. 

2 =2 

0.05271 

0.04726 

0.04150 

0.04369 

1=3» 

j =3 

0.11650 

0.13095 

0.14270 

0.14600 


1 =1, 

2 =2 

-0.08224 

-0.08429 

-0.08610 

-0.09209 


i=l» 

j=l 

0.15980 

0.17530 

0.18870 

0.20080 


i =2» 

j =2 

0.03265 

0.02558 

0.01842 

0.01689 

HS64B 

i=3» 

2 =3 

0.07240 

0.07840 

0.08198 

0.07798 


i 

2 =2 

-0.04889 

-0.04524 

-0.04150 

-0.04038 


i “1 y 

j =1 

0.49410 

0.64560 

0.80250 

0.96740 

HS64G 

i =2» 

j =2 

0.06029 

0.05133 

0.04100 

0.04025 

i =3f 

j =3 

0.17620 

0. 19430 

0.21100 

0.21700 


i=l» 

j =2 

-0. 12570 

-0.13374 

-0.13829 

-0.14679 


1 =1 » 

j =1 

0.07135 

0.07475 

0.07931 

0.08491 

HS64J 







i=2» 

2 =2 

0.02326 

0.01816 

0.01428 

0.01361 


i =3» 

j =3 

0.04360 

0.04696 

0.04753 

0.04428 


i=l» 

j =2 

-0.02565 

-0.02334 

-0.02144 

-0.02099 


i =1 » 

j =1 

0.04595 

0.04829 

0.05160 

0.05602 

HS64M 

i =2f 

j =2 

0.01663 

0.01380 

0.01211 

0.01187 


i =3i 

j =3 

0.03022 

0.03129 

0.03168 

0.03097 


i =1 » 

2 =2 

-0.01708 

-0.01592 

-0.01512 

-0.01535 


i =1» 

j =1 

0.11160 

0.11920 

0.12920 

0.14100 


i =2f 

J =2 

0.04978 

0.04680 

0.04717 

0.05120 

HS640 







i =3» 

J =3 

0.07717 

0.08100 

0.08603 

0.09166 


i =:l» 

i =2 

-0.04044 

-0.04087 

-0.04278 

-0.04612 
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TABLE C-3 


Dissipation Tensor 
Ectuation (4-7) 


Simulation 

EQUATION 

St=4 

St=5 

St=6 

St = 7 


i = l f 


0,43040 

0.49920 

0.55380 

0.60440 


1=2 f 

j=2 

0,05560 

0.05220 

0.05204 

0.05352 

HS64A 







i = 3 » 

j=3 

0,18396 

0.18280 

0.18514 

0.18714 


1=1 f 

j=2 

-0,10240 

-0.10040 

-0.09895 

-0.09735 


1=1 » 

j = l 

0,40380 

0,41180 

0.42000 

0,43420 


i=2» 

3=2 

0,04466 

0.03322 

0.02662 

0.02350 

HS64B 








1=3 » 

j=3 

0.19286 

0.17758 

0.16624 

0.16634 


1=1 » 

j=2 

-0.09403 

-0.07925 

-0.06819 

-0.05889 


i=l t 

j = l 

0,90380 

1.17700 

1.44540 

1.73120 


i=2, 

j=2 

0,05318 

0.04346 

0.03804 

0.03522 

HS64G 






i=3, 

3 = 3 

0.35980 

0.35300 

0.35480 

0.35740 


i=l » 

3=2 

-0.17890 

-0.18600 

-0.19080 

-0.19480 


1=1 f 

3 = 1 

0.07158 

0.06688 

0.06488 

0.06400 


1=2, 

3 = 2 

0,01742 

0.01299 

0.01071 

0.00953 

HS64J 






i=3, 

3=3 

0.05168 

0.04582 

0,04146 

0.03740 


i=l. 

3 = 2 

-0.01863 

-0.01485 

-0.01229 

-0.01054 


1=1 » 

3 = 1 

0.02746 

0.02486 

0,02398 

0.02374 

HS64M 

1=2, 

3 = 2 

0,00940 

0.00680 

0.00564 

0.00508 

1=3, 

3=3 

0,02304 

0.01943 

0.01726 

0.01545 


i=l , 

3=2 

-0.00725 

-0.00568 

-0,00465 

-0.00401 


1=1, 

3=1 

0.06126 

0.05868 

0.05822 

0.05858 

HS640 

1=2, 

J=2 

0.02830 

0.02500 

0.02374 

0.02418 

1=3, 

j=3 

0.05148 

0.04666 

0.04562 

0.04522 


1=1 , 

3=2 

-0.01565 

-0.01432 

-0.01387 

-0.01310 
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TABLE C-4 

Pressu re-strain » Fast Terms 
Eouation (4-46) 


Simulation 

Eouation 

St=4 

St=5 

St=6 

St=7 


i =lf j =1 

-0.17800 

-0.15840 

-0.12280 

-0.12870 

HS64A 

i =2» j =2 

-0.02782 

-0.04524 

-0.03911 

-0.04536 


i =3» j =3 

0.20590 

0.20370 

0.16200 

0.13340 


1 =1 , j =2 

0.27786 

0.24533 

0.16045 

0.13980 


1 =1» j =1 

-0.12380 

-0.10185 

-0.06916 

-0.06710 

HS64B 

1 =2, j =2 

-0.01748 

-0.03175 

-0.02452 

0.00236 

1 =3 » j =3 

0.14134 

0.13360 

0.09368 

0.064*74 


1 =lf J =2 

0.20186 

0.16937 

0.09913 

0.07299 


i J =1 

-1.08190 

-1.00550 

-0.76612 

-0.81830 

HS64G 

i ==2» j =2 

-0.16999 

-0.21173 

-0.19234 

0.06255 

i =3f j =3 

1.25190 

1.21730 

0.95850 

0.75580 


1 =lr j =2 

1.77020 

1.61540 

1 . 12534 

1.00770 


1 =1 » j =1 

-0.02569 

-0.02006 

-0.01378 

-0.01316 

HS64J 

1 =2» j =2 

-0.00394 

-0.00732 

-0.00458 

-0.00070 

1 =3» j =3 

0.02962 

0.02736 

0.01835 

0.01386 


1 j =2 

0.04423 

0.03468 

0.02242 

0.01839 


1 j =1 

-0.01037 

-0.00865 

-0.00671 

-0.00599 

HS64M 

1 =2, j =2 

-0.00147 

-0.00205 

-0.00171 

-0.00101 

1 =3, j =3 

0.01184 

0.01070 

0.00842 

0.00700 


1 ==1» j =2 

0.01707 

0.01482 

0.01130 

0.00960 


1 =1, 2 =1 

-0.03178 

-0.03053 

-0.03006 

-0.03238 

HS640 

1 ==2, j =2 

-0.00386 

-0.00462 

-0.00295 

-0.00332 


i ==3» j =3 

0.03565 

0.03515 

0.03302 

0.03570 


1 =1, j =2 

0.06338 

0.05923 

0.05570 

0.05798 
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TABLE C-5 

Pressure-Strain? Rotta Terms 
Eauation (4-53) 


Simulation 

Eouation 

St=4 

St=5 

St = 6 

St = 7 


i ==1 f 


-0.04969 

-0.06133 

-0.07108 

-0.07963 

HS64A 

i=2r 

j =2 

0.03892 

0.04890 

0.05630 

0.06261 

i=3» 

j=3 

0.01073 

0.01238 

0.01475 

0.01699 


i=l» 

j=2 

0.06990 

0.08400 

0.09830 

0.11640 


i =1 f 

j =1 

-0.01629 

-0,01727 

-0.01722 

-0.01736 

HS64B 

i =2» 

j=2 

0.01338 

0.01506 

0.01543 

0.01586 

i=3. 

j =3 

0.00291 

0.00221 

0.00179 

0.00150 

• 

i=l, 

j=2 

0.02320 

0.02340 

0.02240 

0.02330 


i =1» 

j =1 

-0.03234 

-0.04556 

-0.05800 

-0.06822 

HS64G 

i =2» 

j =2 

0.02912 

0.04139 

0.05148 

0.05567 

i =3, 

j =3 

0.00322 

0.00417 

0.00651 

0.01255 


i=l. 

3=2 

0.05210 

0.06800 

0.07730 

P. 08730 


i =1 f 

j =1 

-0.00698 

-0.00715 

-0.00691 

-0.00676 

HS64J 

i =2, 

j =2 

0,00627 

0.00688 

0.00679 

0.00644 

i =3, 

j =3 

0.00072 

0.00068 

0.00013 

0.00013 


i =1 f 

j =2 

0.01002 

0.00975 

0.00992 

0., 01096 


i =1 » 

j =1 

-0.00348 

-0.00369 

-0.00378 

-0.00396 

HS64M 

i =2f 

3 =2 

0.00320 

0.00343 

0.00373 

0.00377 

i =3f 

3 =3 

0.00027 

0.00025 

0.00005 

0.00019 


1 = 1 » 

3 =2 

0.00587 

0.00568 

0,00555 

0.00620 


i =1» 

3 =1 

-0.01549 

-0.01716 

-0.02000 

-0.02358 

HS640 

1 =2» 

3 =2 

0.01393 

0.01606 

0.01869 

0.02144 

i =3» 

j =3 

0.00155 

0.0011 

0.00130 

0.00214 


i =1 f 

3 =2 

0.03340 

0.03152 

0.03267 

0.03784 
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TABLE C-6 

Scalar Fluctuation Invariance 
2 

< 0 > 


Simulation 

St = 4 

St=5 

St = 6 

St = 7 

HS64A 

0,54163 

0.60711 

0.65946 

0.69973 

HS64B 

0.33816 

0.36453 

0.37910 

0.38215 

HS64G 

0.97368 

1.18150 

1.37179 

1.55312 

HS64J 

0.17704 

0.18678 

0.18985 

0.18687 

HS64M 

0.10957 

0.11203 

0.11355 

0.11522 

HS640 

0.27862 

0.28854 

0.30051 

0,31448 

HS64C 

0.54059 

0.58781 

0.61802 

0.62865 

HS64D 

0.74447 

0.86420 

0.97217 

1.06741 

HS64E 

0.25490 

0.28871 

0,31474 

0.32945 

HS64F 

0.12294 

0.13980 

0.14577 

0.13915 

HS64H 

1.03400 

1.26380 

1.47800 

1.69150 

HS64I 

0.76096 

0.90115 

1.02008 

1.12416 

HS64N 

0.03820 

0.04333 

0.04716 

0.04942 

HS64P 

0.11670 

0.13284 

0.14692 

0.16191 
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TABLE C-7 

Turbulent Heat Flux 


Simulation 

I 

St = 4 

St=5 

St=6 

St = 7 

HS64A 

1 

0.3641 

0.4099 

0.4482 

0.4512 


2 

0.1688 

0.1615 

0.1490 

0.1457 


1 

0.2138 

0.2191 

0.2166 

0.2036 

HS64B 

2 

0.1111 

0.0997 

0.0836 

0.0728 

HS64G 

1 

0.7797 

0.9314 

1.0448 

1 .1412 

2 

0.2674 

0.2698 

0.2604 

0.2588 

HS64J 

1 

0.1359 

0.1376 

0.1338 

0.1265 

2 

0.0625 

0.0539 

0.0433 

0.0371 

HS64M 

1 

0.0775 

0.0797 

0.0832 

0.0823 

2 

0.0397 

0.0341 

0.0292 

0.0273 

HS640 

1 

0.2028 

0.2113 

0.2238 

0.2253 

2 

0.0946 

0.0905 

0.0875 

0.0897 

HS64C 

1 

0.3274 

0.3388 

0,3318 

0.3278 

2 

0.1254 

0.1129 

0.0958 

0.0848 

HS64D 

1 

0.4534 

0.5049 

0.5390 

0.5818 

2 

0.1822 

0.1752 

0.1627 

0.1587 

HS64E 

1 

0.2455 

0.2954 

0,3438 

0.3943 

2 

0.1303 

0.1246 

0.1138 

0.1117 

HS64F 

1 

0.1200 

0.1357 

0.1473 

0.1585 

2 

0.0766 

0.0697 

0.0565 

0.0487 

HS64H 

1 

0.7394 

1.0075 

1.2785 

1.5583 

2 

0.2726 

0.2760 

0.2676 

0.2673 

HS64I 

1 

0.6336 

0.8381 

1.0363 

1.2338 

2 

0.2463 

0.2451 

0.2323 

0.2275 

HS64N 

1 

0.0352 

0.0408 

0.0454 

0.0506 

2 

0.0281 

0.0258 

0.0232 

0.0223 

HS64P 

1 

0.0887 

0.1031 

0,1167 

0.1324 

2 

0.0764 

0.0750 

0.0741 

0.0783 
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TABLE C-8 

Dissipation in the Heat Flux Eauations 
Eouation (4-14) 


Simulations 

i 

St=4 

IH9i 

St=6 

St=7 


1 

0,4079 

0,4862 

0.5546 

0.6103 

HS64A 

2 

-0,1290 

-0,1102 

-0,0983 

-0,0891 


1 

0,3906 

0.4213 

0,4421 

0.4686 

HS(54B 

2 

-0,1274 

-0,0965 

-0.0768 

-0.0611 


1 

0.9269 

1,2450 

1.5495 

1.8530 

HS64G 

2 

-0,2409 

-0.2307 

-0.2229 

-0.2192 


1 

0,0619 

0,0732 

0,0661 

0,0670 

HSI34J 

2 

-0,0272 

-0,0201 

-0.0116 

-0.0082 

HS64M 

1 

0.0214 

0.0221 

0.0227 

0.0234 

2 

-0,0117 

-0.0069 

-0.0040 

-0,0028 


1 

0.0386 

0,0406 

0.0431 

0,0448 

HS640 

2 

-0.0294 

-0.0241 

-0.0201 

-0,0172 

HS64C 

1 

0.3207 

0.3566 

0.3829 

0.3965 


2 

-0,0996 

-0.0770 

-0.0622 

-0.0512 

HS64D 

1 

0.3223 

0,4053 

0.4426 

0.4967 


2 

-0.1061 

-0,0877 

-0,0763 

-0.0646 


1 

0.5024 

0.5364 

0,5924 

0,5717 

HS64E 





2 

-0,1808 

-0,1458 

-0.1322 

-0.1116 


1 

0,4374 

0.4493 

0,4537 

0.4498 

HS64F 

2 

-0,1807 

-0,1380 

-0.1111 

-0.0925 

HS64H 

1 

0,5809 

0,7876 

0.9889 

1.1932 

2 

-0.1500 

-0.1450 

-0.1415 

-0.1405 



2.2700 

2,9279 

3.5175 

4.0632 

HS64I 

2 

-0,6073 

-0,5607 

-0,5225 

-0.4965 


1 

0.0264 

0,0265 

0.0269 

0.0280 

HS64N 

2 

-0,0192 

-0,0137 

-0,0106 

-0,0092 


1 

0,0645 

0.0669 

0,0700 

0.0737 

HS64P 

2 

-0,0571 

-0,0489 

-0.0431 

-0.0407 
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TABLE C-9 

Pressure-Scalar Covariance) Fast Term 
Eouation (4-60) 


Simulation 

i 

St=4 

St=5 

St = 6 

St = 7 

HS64A 

1 

-0«5182 

-0.4350 

-0.2959 

-0.2565 

2 

0.1842 

0.1986 

0.1740 

0.0768 

HS64B 

1 

-0.3897 

-0.3136 

-0.1923 

-0.1448 

2 

0.1414 

0.1592 

0.1369 

0.0597 

HS64G 

1 

-3.2060 

-2.8110 

-1.9854 

-1.7207 

2 

1.6390 

1.6950 

1.5810 

0.9050 

HS64J 

1 

-0.0832 

-0.0630 

-0.0400 

-0.0309 

2 

0.0227 

0.0274 

0.0207 

0.0096 


1 

-0.0310 

-0.0254 

-0.0194 

-0.0164 

HS64M 

2 

0.0078 

0.0076 

0.0063 

0.0042 


1 

-0.1044 

-0.0981 

-0.0932 

-0.0932 

HS640 

2 

0.0188 

0.0179 

0.0151 

0.0145 

HS64C 

1 

-0.3353 

-0.2312 

-0.1194 

-0.0816 

2 

0.1297 

0.1209 

0.0871 

0.0333 

HS64D 

1 

-0.5344 

-0.4475 

-0.3061 

-0.2662 


0.1938 

0.2040 

0.1759 

0.0751 

HS64E 

1 

-0.4605 

-0.3923 

-0.2641 

-0.2283 

2 

0.1501 

0.1763 

0.1614 

0.0741 

HS64F 

1 

-0.3174 

-0.2618 

-0.1565 

-0.1165 

2 

0.0966 

0.1230 

0.1085 

0.0442 


1 

-3.2840 

-2.8303 

-2.0009 

-1.7348 

HS64H 

2 

1.6716 

1.7254 

1.6090 

0.9296 

HS64I 

1 

-3.1525 

-2.7229 

-1.9157 

-1.6583 

2 

1.4977 

1.5639 

1.4650 

0.8042 

HS64N 

1 

-0.0234 

-0.0205 

-0.0161 

-0.0140 

2 

0.0055 

0.0063 

0.0059 

0.0042 


1 

-0.0865 

-0.0831 

-0.0786 

-0.0796 

HS64P 

2 

0.0155 

0.0168 

0.0155 

0.0144 
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TABLE C-1.0 


Pressure-Scalar Covariance} Rotta Term 
Eouation (4-74) 


Simulation 


HS64A 


HS64B 


HS64G 


HS64J 


HS64M 


HS640 


HS64C 


HS64D 


HS64E 


HS64F 


HS64H 


HS64I 


HS64N 


HS64P 












St = 4 

St=5 

St=6 

St=7 

-0.1209 

-0.1231 

-0.1258 

-0.1314 

0.1816 

0.1936 

0.1964 

0.2082 

-0.0444 

-0.0379 

-0.0324 

-0.0300 

0.0687 

0.0635 

0.0542 

0.0522 

-0.0812 

-0.0865 

-0.0946 

-0.1062 

0.1502 

0.1664 

0.1573 

0.1694 

-0.0222 

-0.0193 

-0.0169 ► 

-0.0144 

0.0296 

0.0281 

0.0269 

0.0268 

-0.0123 

-0.0115 

-0.0096 

-0.0085 

0.0157 

0.0158 

0.0146 

0.0136 

-0.0601 

-0.0612 

-0.0648 

-0.0692 

0.0711 

0.0696 

0.0732 

0.0772 

-0.0467 

-0.0352 

-0.0254 

-0.0205 

0.0731 

0.0590 

0.0437 

0.0343 

-0.1417 

-0.1473 

-0.1538 

-0.1625 

0.2125 

0.2303 

0.2381 

0.2678 

-0.0678 

-0.0664 

-0.0673 

-0.0707 

0.0997 

0.1010 

0.0981 

0.1009 

-0.0176 

-0.0141 

-0.0124 

-0.0110 

0.0257 

0.0217 

0.0175 

0.0142 

-0.0839 

-0.0902 

-0.0988 

-0.1114 

0.1547 

0.1727 

0.1635 

0.1767 

-0.0693 

-0.0710 

-0.0771 

-0.0859 

0.1298 

0.1383 

0.1297 

0.1370 

-0.0055 

-0.0054 

-0.0047 

-0.0043 

0.0061 

0.0065 

0.0063 

0.0062 

-0.0365 

-0.0378 

-0.0450 

-0.0458 

0.0426 

0.0425 

0.0465 

0.0500 
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